宏基因组reads筛选:去除宿主序列

原创
2019/12/23 11:36
阅读数 929

基于环境的复杂性与研究对象的不同,宏基因组数据在组装之前常需要过滤掉一些序列以防干扰研究。例如要研究动植物组织或肠道的微生物组,往往需要去除宿主的DNA序列。假如研究的是人类肠道微生物的宏基因组,需要去除属于人基因组的序列。具体方法为将质控后的序列和人类基因组序列进行比对,将比对上的序列去除。

宏基因组分析Pipeline

测序数据的解析: Fastq与FastQC
测序数据的质控: Trimmomatic!
宏基因组reads筛选: 去除宿主序列
测序数据的组装:常用软件工具
更新中……
短序列有参比对常用的软件有 BWA Bowtie BBMap 等。下面以 Bowtie 2 为例。 Bowtie 2 http://bowtie-bio.sourceforge.net/bowtie2/index.shtml )是一个快速的、节约内存的序列比对工具,用于将测序的 reads 比对到长的参考序列。首先需要下载参考基因组,这里以人类为例,可以去 NCBI 下载最新版本的人类基因组序列( https://www.ncbi.nlm.nih.gov/projects/genome/guide/human/index.shtml ),其下载的为 fasta 格式(压缩文件),如下所示:

染色体两端为端粒重复序列所以用 N 标记,接下来解压文件然后使用 bowtie2-build 来构建新的 index ,如下所示:
gzip -d GRCh38_latest_genomic.fna.gzbowtie2-build --threads 20 GRCh38_latest_genomic.fna human_genome
运行结束后,生成 6 个文件,如下所示:

Bowtie 已经为很多常见物种预先构建了 index (见官网主页 Pre-built indexes 栏),不想自己构建可以在其 FTP 中下载。 bowtie2 进行比对的用法如下所示:
bowtie2 [options] -x <bt2-idx> { -1 <m1> -2 <m2> | -U <r>} [-S <hit>]<文件>-x <bt2-idx>参考基因组(reference genome)通过bowtie2-build指令构建的Index文件-1 <m1>双末端测序中第一个fastq文件,可以写多个文库但是必须用逗号隔开,但文件m1与文件m2必须一一对应,测序文件中的Reads的长度可以不同。-2 <m2>双末端测序对应的第二个fastq文件,与文件m1对应-U <r>与前面的文件1,文件2为或的关系,此处的文件是非双末端比对文件。例如lane1.fq,lane2.fq,lane3.fq,lane4.fq。可以是多个文件,但是必须用逗号隔开。-S <hit>指定输出文件,后缀是sam的格式的文件,默认标准输出[options]:-qReads(用<m1><m2><s>指定)是FASTQ格式的文件,默认即FASTQ。--qseqReads(用<m1><m2><s>指定)是QSEQ格式的文件。-fReads(用<m1><m2><s>指定)是FASTA文件。-rReads(用<m1><m2><s>指定),每行代表一个输入序列,没有任何其他信息(无read名称,无qualities)。-c后面直接是比对的reads序列(而不是文件),即reads序列在命令行上给出。-s/--skip <int><int>中是数字,input的reads跳过前<int>个reads或read pairs-u/--qupto <int>比对前<int>个reads或read pairs,然后停止。-5/--trim5 <int>剪掉5'(左)端<int>长度的碱基,再用于比对(默认值0)-3/--trim3 <int>剪掉3'(右)端<int>长度的碱基,再用于比对(默认值0)--phred33输入的序列质量数据为Phred33体系(默认为phred33)--phred64输入的序列质量数据为Phred64体系-p程序运行所用核数
接下来进行比对,如下所示:
bowtie2 -p 20 -x human_genome -1 clean_1.fq -2 clean_2.fq -S meta.sam
运行结束后生成 sam 格式的结果文件,用来存储 reads 到参考序列的比对信息(即告诉你这个 reads 在染色体上的位置等)。接下来使用 samtools 工具将 sam 文件转化为 bam 文件,其中 -b 指定输出文件为 bam -S 指定输入文件为 sam
samtools view -bS --threads 20 meta.sam > meta.bam
bam 文件为 sam 文件的压缩版本,是二进制文件,占用空间更小,可以用 samtools 工具实现 sam bam 文件之间的转化。接下来对 bam 文件按照比对的位置坐标对 reads 进行排序:
samtools sort meta.bam -o meta.reads.sorted.bam --threads 20

继续将bam文件转换为bed文件:

bamToBed -i meta.reads.sorted.bam > meta.reads.sorted.bed

bed文件中包含了全部比对到宿主基因组的序列信息,根据序列信息,将原始数据中包含有宿主基因组的序列去除

其中第一列为参考基因组染色体或 scaffold 名称,第二列与第三列为 read 在该染色体或 scafflold 比对的起始与终止位置,第四列为比对上的 read 名称(末尾 /1 表示第一端 /2 为第二端,一般只要有一端比对上,双端的序列均要去除),第五列为得分 score ,第六列为比对到参考序列链的方向。
最后,需要大家自行写脚本将这些reads从测序数据中去除。
END

本文分享自微信公众号 - 微生态与微进化(MicroEcoEvo)。
如有侵权,请联系 support@oschina.cn 删除。
本文参与“OSC源创计划”,欢迎正在阅读的你也加入,一起分享。

展开阅读全文
打赏
0
0 收藏
分享
加载中
更多评论
打赏
0 评论
0 收藏
0
分享
返回顶部
顶部