paired-end reads的拼接

paired-end reads的拼接

发表于2012 年 8 月 13 日

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

参考：

• http://www.ebi.ac.uk/~zerbino/velvet/
• http://www.bigsa.org.au/node/92
• http://seqanswers.com/wiki/ABySS#Paired-end_assembly