分享

更多

   

FASTA序列文件处理一网打尽

2017-08-09  海陵岛的鱼

推荐两个地方:

地方一都是小脚本,但实用,大伙也可以自己练习写。

地方二成熟软件SeqKit,也很实用。

一、小脚本

大家可以在这里下载以下脚本:

https://github.com/jorvis/biocode/tree/master/fasta


各脚本作用信息如下:


|-- append_to_fasta_header.py

每个序列ID添加后缀

|-- check_for_embedded_fasta_headers.py

检查此种错误类型fasta序列

>gi1006569 DnaA [Vibrio harveyi]

MSSSLWLQCLQQLQEELPATEFSMWVRPLQAEVLHAC>gi409247 DNA repair protein

MVSLTFKNFKKEKVPLDLEPSNTILETKTKLAQSISCEESQIKLIYSGKVLQDSKTVSECGLKDGDQVVF

|-- compare_two_fastas.pl

比较两个文件相同序列与不同序列数目

|-- convert_fasta_contigs_to_gff3.py

contig序列转换为GFF格式

|-- convert_fastq_to_fasta.py

FASTQ转换为FASTA

|-- create_fasta_pseudomolecules.pl

从遗传图或者Hi-C确定的map信息文件中获得假染色体序列(contig连成染色体序列)

|-- extract_fasta_regions.py

提取特定区域的序列

|-- fasta_base_content.py

碱基含量统计

|-- fasta_size_distribution.pl

GC含量和长度分布统计

|-- fasta_size_distribution_plot.py

序列长度分布作图

|-- fasta_size_report.pl

序列长度信息统计

|-- filter_fasta_by_header_regex.py

使用正则表达式获取FASTA序列

|-- filter_fasta_by_ids.pl

根据ID提取序列

|-- filter_fasta_by_size.pl

根据长度过滤序列

|-- merge_fasta_files_and_uniquify_ids.py

你有多个序列ID相同的文件,想合并一起,保证序列ID唯一性

 >CL1Contig1

 >CL1Contig2

|-- merge_masked_fasta_files.py

将不同软件冰壁重复序列结果合并

|-- prepend_to_fasta_header.py

你有多个序列ID相同的文件,想合并一起,在序列前面添加前缀,保证序列ID唯一性

|-- reformat_fasta_residue_lengths.py

按照一行固定字符数目重新显示FASTA序列

|-- remove_duplicate_sequences.py

移除掉重复的序列

|-- remove_empty_sequences.pl

移除掉文件中空的序列

|-- reorient_fasta.pl

从某一位置截取序列到另一个文件中

|-- reorient_sequences_by_id.py

将指定序列反向或者反向互补

|-- split_fasta_into_even_files.py

分割成指定数目的文件

|-- split_interleaved_sequence_file.py

FASTQ/FASTA序列中成对的分成R1,R2两个文件,单独非成对的生成一个文件

|-- subsample_fasta.py

随机抽取指定数目的FASTA序列

|-- validate_fasta.py

主要检查以下内容:

- ERROR: Entry with 0-length sequence

- ERROR: No > symbols embedded in sequence residues (suggests merged records)

- ERROR: Nonzero-length identifier after > symbol until first whitespace

- ERROR: Multiple records within an individual file with same identifier

- ERROR: Incorrect record count (if --expected_record_count is passed)

- ERROR: Make sure the file ends with a blank line (else commonly used tools like cat will cause errors)

- WARNING: Residue lines > 60bp

- WARNING: Comment lines with # (most parsers will fail to handle this)

- WARNING: Homopolymers of length > --homopolymer_limit which are not Ns

- WARNING: Internal stops (*) detected in protein FASTA (if --check_internal_stops is passed)

`-- write_fasta_from_gff.pl

GFF提取CDS和蛋白序列

 

二、SeqKit


       处理FASTA和FASTQ的神器,window\linux系统版本都有.对于没有编程基础的小伙伴们,我们照样可以轻松操作序列文件。

该软件功能强大,小编只罗列部分模块功能,更详细功能参见软件网站:

http://bioinf.shenwei.me/seqkit/usage/。


一、序列操作。


  seqkit seq [flags] file

参数:

  -p, --complement           取互补序列

      --dna2rna                   DNA to RNA

  -l, --lower-case                将序列以小写字母形式输出

  -g, --remove-gaps           移除组装序列中的gap

  -r, --reverse                    取反向序列

  --rna2dna                     RNA to DNA

  -u, --upper-case                将序列以大写字母形式输出

-w, --line-width int               以每行指定长度输出序列 (0 for no wrap) (default 60)

举例:

seqkit seq test.fa -w 0#将此文件fasta序列转换成一行输出

seqkit seq -w 100  test.fa#将此文件fasta序列转换成100个碱基一行输出

seqkit seq --dna2rna test.fa#将此文件fasta序列dna转换成rna

seqkit seq -w 100 -p -r test.fa#将此文件fasta序列反向互补输出,每行100碱基



二、Fasta/q之间及与tab格式互换



1、FASTQ转换成FASTA: seqkit fq2fa 

举例:

seqkit fq2fa reads_1.fq -o reads_1.fa

2、FASTA/FASTQ转换成tab格式。seqkit fx2tab

举例:

seqkit fx2tab test.fa>test.fa.tab.fa

seqkit fx2tab test.fq>test.fq.tab.fq

tab格式:ID  sequence




三、序列信息统计



1、序列碱基含量及序列长度信息统计

  seqkit fx2tab [flags]

参数:

  -B, --base-content value   要输出的碱基含量e.g. -B AT -B N

  -g, --gc                   print GC content

  -l, --length               print sequence length

  -n, --name                only print names

-i, --only-id              print ID instead of full head

举例:

seqkit fx2tab  -l -g -n -i -H  test.fa 

输出结果:

#name seq     qual     length     GC

gene1                      30                 40.00

2、序列长度分布统计

Usage:

seqkit stat [flags]

举例:

seqkit stat test.fa

输出结果    

file     format  type  num_seqs  sum_len  min_len  avg_len  max_len

test.fa  FASTA   DNA          1       30       30       30       30




四、根据ID或特定的motif筛选提取序列




seqkit grep [flags]

参数:

  -n, --by-name               匹配整个序列的名字,包含description部分,而不是序列id。  

  -s, --by-seq                 匹配序列

  -d, --degenerate            pattern/motif 包含简并碱基

  -i, --ignore-case            忽略大小写

  -v, --invert-match          输出不匹配此模式的内容

  -p,            匹配模式,支持连续写多个模式,匹配任一模式即输出。如-p ^ATG -p TAA$。注意该功能仅能正向匹配,不能实现对互补链匹配。

  -f, --pattern-file string   支持匹配模式写到一个文件中,如要提取的序列ID。

  -R, --region string         匹配位置选择。e.g 1:12 for first 12 bases, -12:-1 for last 12 bases

  -r, --use-regexp            使用正则表达式,必须加入此参数,如^匹配首端。同-p联合使用。

示例:

seqkit grep -s -r -i -p ^atg cds.fa#选取有起始密码子的序列

seqkit grep -f list test.fa > new.fa#根据ID提取序列

seqkit grep -s -d -i -p TTSAA#简并碱基使用。S 代表C or G.

seqkit grep -s -R 1:30 -i -r -p GCTGG##匹配限定到某区域



五、motif定位



对grep的拓展,可以正反链同时匹配,输出匹配的位置。

seqkit locate [flags]

参数

  -d, --degenerate                pattern/motif contains degenerate base

  -i, --ignore-case                ignore case

  -P, --only-positive-strand      only search at positive strand

  -p, --pattern value             search pattern/motif 

  -f, --pattern-file string       pattern/motif file (FASTA format)

举例

seqkit locate -i -d -p AUGGACUN test.fa

输出结果

seqID         patternName   pattern    strand   start   end   matched

cel-mir-58a   AUGGACUN      AUGGACUN   +        81      88    AUGGACUG

ath-MIR163    AUGGACUN      AUGGACUN   -        122     129   AUGGACUC




六、多个序列文件比较寻找相同的序列或者ID相同的序列




 seqkit common [flags]

参数:

-n, --by-name       匹配整个序列的名字,包含description部分,而不是序列id 

-s, --by-seq         match by sequence

-i, --ignore-case   ignore case

-m, --md5           use MD5  reduce memory usage 

举例:

1、By ID (default,>后面,空格之前的名字)输出ID名字相同的。

seqkit common test1.fa test2.fa -o common.fasta

2、By full name(整个序列的名字,包含description部分)。输出序列名字相同的。

seqkit common test1.fa test2.fa  -n -o common.fasta

3、输出要比较的文件中序列相同的序列

seqkit common test1.fa test2.fa -s -i -o common.fasta

4、输出要比较的文件中序列相同的序列 (for large sequences)

seqkit common test1.fa test2.fa -s -i -o common.fasta --md5



七、提取部分序列



如随机抽取10000条FASTQ序列做NT污染评估。同时他也可以对FASTA序列提取

seqkit sample [flags]

参数:

  -n, --number int         sample by number (result may not exactly match)

  -p, --proportion float   sample by proportion

  -s, --rand-seed int      rand seed for shuffle (default 11)

  -2, --two-pass           2-pass mode lower memory

 举例:随机抽取序列

 seqkit sample -n 10000  -s 11  test1_1.fq -o sample.fq

seqkit sample -p 0.1 -s 11  test1_1.fq -o sample.fq



八、排序输出命令



seqkit sort [flags]

参数:

  -l, --by-length               按照序列长度排序

  -n, --by-name                 by full name 

  -s, --by-seq                  按照序列排序

  -i, --ignore-case            按序列排序时忽略大小写

  -r, --reverse                 反向排序

  -2, --two-pass                对于FASTA序列排序可以减少内存

举例:

seqkit sort -l test.fa




九、文件切割




seqkit split [flags]

参数:

  -i, --by-id              split squences according to sequence ID

  -p, --by-part int        将一个文件分割成N 份

  -s, --by-size int        将一个文件按照N 条序列一个文件进行分割

  -O, --out-dir string     output directory (default value is infile.split)

  -2, --two-pass           two-pass mode  to lower memory usage(only FAST)

举例:

seqkit split hairpin.fa.gz -p 4



    本站是提供个人知识管理的网络存储空间,所有内容均由用户发布,不代表本站观点。如发现有害或侵权内容,请点击这里 或 拨打24小时举报电话:4000070609 与我们联系。

    猜你喜欢

    0条评论

    发表

    类似文章 更多
    喜欢该文的人也喜欢 更多