paired-end reads的拼接

Velvet中paired-end reads的拼接

文件格式

要将两头测序(paired-end)的reads放到同一个文件当中,fastq格式,必须成对的依次放置reads [interleaved],velvet是成对读取的,另外Velvet假设来自两头read是反向互补的,如果不是,需要用反向互补序列来代替第一个read。Fastq格式中paired-end reads的编号相同,但是其有/1或者/2的后缀,通过这种方式来标示paired-end reads。

如果两端测序的reads放在不同的两个文件中,可以使用Velvet提供的perl脚本shuffleSequences fasta.pl进行转换合并,命令格式如下:

> ./shuffleSequences_fasta.pl forward_reads.fa reverse_reads.fa output.fa

低质序列过滤

在拼接前,首要要进行去除低质序列、接头等预处理,比如使用FASTX-Toolkit中的fastq_quality_filter去除低质序列:

fastq_quality_filter  -q 20 –p 100 -i s_1_1_sequence.txt -o s_1_1_sequence.txt_filtered_q20_p100.fastq
fastq_quality_filter  -q 20 –p 100 -i s_1_2_sequence.txt -o s_1_2_sequence.txt_filtered_q20_p100.fastq

这样势必带来一个问题,有些paired-end的前面序列被剔除,有些后面的序列被剔除,paired-end序列无法成对的错落出现,下面需要做的就是必须将单独的reads挑出来,方法有很多,下面是其中一个:

合并到一个文件中

cat s_1_[12]_sequence.txt_filtered_q20_p100.fastq > s_1_filtered_q20_p100.fastq
rm s_1_[12]_sequence.txt_filtered_q20_p100.fastq

使用cdbfasta为Fastq创建索引

cdbfasta -Q s_1_filtered_q20_p100.fastq

导出所有序列编号

cdbyank s_1_filtered_q20_p100.fastq.cidx -l > s_1_filtered_q20_p100.fastq.ids

使用awk根据序列编号的特点,/1或者/2后缀,对于编号进行过滤

#得到完整paired-end reads
awk -v sep="/" '{ if ((sep_i=index($0,sep)) > 0) { name=substr($0,1,sep_i-1); suffix=substr($0,sep_i); } else { name=$0; } if (r[name]) { print name r[name]; print $0; delete r[name]; } else { r[name]=suffix; }}' s_1_filtered_q20_p100.fastq.ids > s_1_filtered_q20_p100.fastq.paired.ids
#得到单独的reads(orphaned reads)
awk -v sep="/" '{ if ((sep_i=index($0,sep)) > 0) { name=substr($0,1,sep_i-1); suffix=substr($0,sep_i); } else { name=$0; } if (r[name]) { delete r[name]; } else { r[name]=suffix; }}END {for (name in r) print name r[name]}' s_1_filtered_q20_p100.fastq.ids > s_1_filtered_q20_p100.fastq.orphans.ids

根据编号,得到相应的Fastq格式的序列文件

cdbyank s_1_filtered_q20_p100.fastq.cidx < s_1_filtered_q20_p100.fastq.paired.ids > s_1_filtered_q20_p100.fastq.paired.fastq
cdbyank s_1_filtered_q20_p100.fastq.cidx < s_1_filtered_q20_p100.fastq.orphans.ids > s_1_filtered_q20_p100.fastq.orphans.fastq

 运行velveth

> ./velveth output_directory/ 21 -fastq -shortPaired s_1_filtered_q20_p100.fastq.paired.fastq -fastq -short s_1_filtered_q20_p100.fastq.orphans.fastq

运行velvetg

 > ./velvetg output_directory/ -ins_length 400 -exp_cov 21.3

使用ABySS拼接

abyss-pe k=25 n=10 in='s_1_filtered_q20_p100.fastq.paired.fastq' se='s_1_filtered_q20_p100.fastq.orphans.fastq' name=my_organism

参考:

paired-end reads的拼接》上有1条评论

  1. 你好,我用fastx_artifacts_filter 来去除低质量的reads,可是报错出Invalid quality score value (char ‘#’ ord 35 quality value -29) on line 4
    ,不能识别’#’的质量值,请问有什么解决的方法吗?

发表评论

电子邮件地址不会被公开。 必填项已用*标注

请启用Javascript,以完成验证!