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