什么是Barcode?
因为测序仪的通量很高,一台测序仪拥有多个数据量可以达到千万级别的文库,并且通过index区分不同的文库。而我们获得16S的数据之后可以发现每个样本的数据量很小(从文件大小就可以看出来),因此实际上我们是将多个样品混合在一起构建成一个测序文库,而不是一个样品单独成为一个测序文库。
那么为了区分在同一个文库中来自不同样品的序列,我们就需要给序列贴上唯一的标识符,这个标识符就是在序列引物外侧的Barcode了。通常barcode的长度为12bp。
原始测序数据有哪些?
通常通过MiSeq仪器测序会获得三个文件。包含'_R1_' 的是forward 序列,包含'_R3_'的是reverse序列,包含'_R2_'的是barcode序列。为了确保无误,可以确认一下含有R2的文件大小是否小于其他两个文件。
如果获得的是来自于公司的文件,通常公司已经进行过分样处理。所谓的分样处理,就是将forward 序列或reverse序列按照样本来源即根据Barcode进行分类并切除引物,最终形成单个样品自己的forward 序列和reverse序列。因此,我们最后获得的可能是每个样本的forward 序列和reverse序列,而没有上述的barcode文件及整合的forward 序列和reverse序列文件。
导入序列文件
我们一般从公司获得都是双端测序已经除去barcode和引物的,按照分样(根据样品分类)完毕的数据。所以第一步我们通过flash软件对R1和R2序列数据进行拼接。
合并序列
#安装flash软件conda install flash#查看使用方法flash -h#拼接序列 2>&1表示将标准错误重新定向到标准输出#大家可以对比有无 >joined/1.log 会发现没有这段代码表示程序运行的结果输出到屏幕,有则生成文件1.log#可以自行修改默认参数具体见flash -hmkdir joined flash rawdata/1_R1_.fastq.gz rawdata/1_R2_.fastq.gz -o 1 -d joined/ >joined/1.log 2>&1 flash rawdata/2_R1_.fastq.gz rawdata/2_R2_.fastq.gz -o 2 -d joined/ >joined/2.log 2>&1
输出结果:
-1.extendedFrags.fastq 拼接后的序列-1.notCombined_1.fastq R1中没有成功拼接的序列-1.notCombined_2.fastq R2中没有成功拼接的序列-1.hist 拼接序列的长度的直方图(Numeric)-1.histogram 拼接序列的长度的直方图(Visual)-1.log 程序运行日志
保留合并后文件和mapping文件
mkdir 1mv joined/1.extendedFrags.fastq map/1.txt 1/ mkdir 2mv joined/2.extendedFrags.fastq map/2.txt 2/
合并后序列质控
#具体参数可以见split_libraries_fastq.py -hsplit_libraries_fastq.py -i 1/1.extendedFrags.fastq -m 1/1.txt -q 30 --barcode_type not-barcoded --sample_id 1 -o 1 --store_demultiplexed_fastq split_libraries_fastq.py -i 2/2.extendedFrags.fastq -m 2/2.txt -q 30 --barcode_type not-barcoded --sample_id 2 -o 2 --store_demultiplexed_fastq
注意:这一步没有-m参数也可以,但是要注意--sample_id这个参数后面与mapping中的要求一致,只能使用“.”和字母数字,否则后续生成的文件会出错。
输出文件:
合并所有的fna文件
#利用cat命令cat */*.fna> seq.fna
导入其他格式的序列文件
导入其他格式的序列文件可以参考qiime1的教程。
作者:jlyq617
链接:https://www.jianshu.com/p/f5f6caa7982a
共同学习,写下你的评论
评论加载中...
作者其他优质文章