|
引子 目前NCBI上的数据已经超过的14个P，公开的测序数据也已近到达5个P，而且目前随着二代测序成本的继续下降，测序数据将在很长一段时间内持续井喷。面对如此大量的数据包含大量的与作物相关土壤相关的数据，就像是面对一大块肥肉一样，这些数据往往没有得到充分的利用，结合目前的在农业上方兴未艾（医学上已近司空见惯）meta分析，相信我们利用这种方式处理大量的测序数据是大数据时代科研的一个重要手段和方向。 所以需要在测序数据分析上达到快速化，脚本化，和流程化。为了这个目的，我的这个系列的推送将围绕这个主题，从NCBI上得到数据开始，依托ubuntu系统，使用qiime和usearch等众多生信分析手段，优化16s分析策略，构建特定功能的一体化分析流程。今天我关注的是：迅速提高序列拼接效率，得到后续分析友好型输入，依托qiime。 我们使用Qiime对二代双端测序数据进行拼接总会遇到这样的尴尬： 目前很多公司交给我们的测序结果都是每个样品双端序列分来的两个fastq文件，在qiime中我们有两种处理方式： 第一就是（命令链接：http:///scripts/join_paired_ends.html）：使用 join_paired_ends.py -f$PWD/forward_reads.fastq -r $PWD/reverse_reads.fastq -o $PWD/fastq-join_joined 这条命令只能指定单个样品的两个正反序列进行拼接，多个样品我们就挨个编命令（这边是我第一次使用qiime所用的方式）。这里限制速度的有两个因素：1.需要逐条编写命令；2.输出结果的拼接好的文件fq命名都是一样的，还需要将这些修改对应的样品名。3。需要将这些单个文件从对应的每个文件夹中提取出来，才方便后续的质控操作。 第二就是（命令链接：http:///scripts/multiple_join_paired_ends.html）使用： multiple_join_paired_ends.py -i input_files-o output_folder -p qiime_parameters.txt 这条命令省去了我们逐条编写每个样品拼接的麻烦，但是输出结果和上一条命令一样，拼接好的fq文件名字都一样，都存在着修改名字和转移这个fq文件的耗费时间的操作。（在过去的一段时间我使用的这个命令） 我的目的是按照自己的分析要求，组合相关的功能，构建完整的16s分析流程，只需要输入数据，静等结果的输出。但是序列拼接这一步需要大量的人工操作，实在是一个障碍，于是今天关注这个障碍： 在早先的usearch11推出来的时候我跟着他的流程跑了一遍，用了几个小时（想起当时第一次学Qiime，两份流程跑了一个月）。上面对双端序列拼接的方式我很感兴趣：这里摘抄下来：这是一个简单的shell的for循环，实现了三个功能，第一是批量的序列拼接，第二就是将样品名添加到序列名形成上，第三就是将所有序列合并，方便后续去冗余。 for Sample in Mock Soil Human Mouse; do usearch11 -fastq_mergepairs${Sample}*_R1.fq -fastqout$Sample.merged.fq -relabel $Sample.; cat $Sample.merged.fq >>all.merged.fq; done 这里我们采用这个方式，依托qiime来进行双端序列拼：首先是建立两个文件夹存储数据；然后我们切换到双端测序结果文件夹中，使用*_1.fastq定位全部测序样品（可能有人的样品后缀是.fq，都一样的，换过来）；第二步，使用{i%_*}来提取样品名称（使用字符串裁剪工具：意思就是去除左边第一个_后的内容），第三步集成join_paired_ends.py -f $PWD/forward_reads.fastq -r$PWD/reverse_reads.fastq -o $PWD/fastq-join_joined这条命令对每个样品进行序列拼接；第四步将所有拼接好的文件转移到一个新的文件中；这里我使用的是-m SeqPrep的方法进行拼接，（据说会好一些，默认不是这个），但是结果是已gz结尾的压缩文件，所以这里我们多加一步，解压缩,并删去原来的压缩文件。 mkdir split mkdir pair for i in *_1.fastq; do echo ${i%_*}; join_paired_ends.py -m SeqPrep -f$PWD/${i%_*}_1.fastq -r $PWD/${i%_*}_2.fastq -o $PWD/split/${i%_*}; mv $PWD/split/${i%_*}/seqprep_assembled.fastq.gz$PWD/pair/${i%_*}.fastq.gz; tar xvzf $PWD/pair/${i%_*}.fastq.gz; rm $PWD/pair/${i%_*}.fastq.gz; done 这样一来我们质控便十分轻松了，省去了大量的人工修改各种东西的操作。 学习永无止尽，分享永无停歇。 |
|
|
