老司机的扩增子分析排坑指南 / 开普饭

这是一篇排坑指南，并不是操作指南，主要是理解思想，出现错误可以找到有限的候选方向去解决。

写在前面：

扩增子分析不知不觉已经做了三年了，什么大坑都见过了，如何修改都经历了，可是今天我做的一组数据竟然完全cover了我三年来做扩增子遇到的大部分坑。一方面这是一个恶心的事情，毕竟我虽然是老手，但是解决问题还是花费了我一天时间，另一方面呢：这个过程记录下来确实可以帮到大家，还有以后的我。话不多说，让我们跟着昨天晚上的思路进行操作。

第一个坑：

两年前我们使用qiime来做分析，我用的虚拟机，这是安装好的环境，但是如今我使用服务器，使用conda构建的虚拟环境，这两者不同，虚拟环境有些功能并没有安装，比如遇到的：join_paired_ends.py 质控并没有两个质控方法模块

for i in *_1.fq; do echo ${i%%_*}; join_paired_ends.py -m SeqPrep -f $PWD/${i%%_*}_1.fq -r $PWD/${i%%_*}_2.fq -o $PWD/resistant_join1/${i%%_*}; mv $PWD/resistant_join1/${i%%_*}/seqprep_assembled.fastq.gz $PWD/resistant_join2/${i%%_*}.fastq.gz; gzip -d $PWD/resistant_join2/${i%%_*}.fastq.gz; done

解决：

因为 SeqPrep质控效果根更好，所以我安装了 SeqPrep模块：

conda install SeqPrep

拿到数据第二个坑：

这批数据是两年前测的，虽然是测序在同一个公司，同一个引物，但是得到数据的我发现这明显是两个平台测定的，以@HISEP和@M000开头的两种数据格式各都有，因此就出现了第一个坑：以@M000开头的数据并不能很好的进行质控：我使用的Qiime质控，他要求严格的序列命名

multiple_split_libraries_fastq.py -i a2_join/ -o a3_resistant_split -p /home/qiime/Desktop/Shared_Folder/mel_split_tmt.txt --demultiplexing_method sampleid_by_file

解决：

于是写了个python脚本进行修改名称：

import osimport reimport datetimedef replace_content(name): with open('a0_raw_data/{name}'.format(name=name), 'r') as f: data = f.read() data1 = re.match('\@\w{3}\d{7}.', data).group() data2 = re.sub(data1,'@HWI-D00477:901:HT3KKBCX2:1:1101:', data) data3 = data2.replace(' ', ':') data4 = re.sub(r':length=251', ' 1:N:0:TATGGCAC', data3)

with open('a0_raw_data/{name}'.format(name=name), 'w') as f: f.write(data4)

global data2names = os.listdir('a0_raw_data')starttime = datetime.datetime.now()for name in names: print(name) replace_content(name)endtime = datetime.datetime.now()duringtime = endtime - starttimeprint(duringtime)

第三个坑：

之前我们测定的V3-V4区域，但是使用的确实PE250，这将表明不可用DADA2，所以你们也就别尝试了，一般而言前端序列质量优于后端，要么仅仅使用前端序列进行分析，要么你就需要免面对拼接率只有一半甚至更少的情况，其实无论是DADA2，常规基于聚类的扩增子分析的拼接率也不高，最多一般，已过后端序列质量差的情况，所以我认为这么长的片段以后就不要测了，要么换平台，要么换引物了。

第四个坑：

可用序列过少问题，这是两年前的测序结果，但是我通过upase聚类得到的可用序列数量不到40%。这种数据只有在NCBI中下载早先的测序数据才会碰到，我觉得很器官，然后我又尝试了参考聚类方法，unoise3聚类，都是同样的情况，因此我往回追溯，质控问题，拼接问题？最终还是发现是拼接问题，拼接成功率不到40%，所以这就没办法了，只能舍弃后端序列。

第五个坑：

因此我想一条命令提取前端序列：

mv *1.fq.gz ./raw/

事实证明我太天真了，有些序列虽然被公司标记为1，但却不是前端，我的亲娘哎，这只能让我一个样本一个样本去NCBI比对，看看序列是哪个区间，本次我分析的数据量为104个样本，共208个前后端数据。我足足比对了好几个小时，发现有三分之一的样本标记是错误的。这个过程让我怀疑了人生。

写在后面

零零散散的一套下来，一晚上加一早上就过去了，没想到吧，这是命运带我的的礼物，让我复习知识呀。

老司机的扩增子分析排坑指南

相关推荐