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