宏基因组实战7. bwa序列比对, samtools查看, bedtools丰度统计
前情提要
如果您在学习本教程中存在困难，可能因为缺少背景知识，建议先阅读本系统前期文章
测试数据
刘博士帮助把测试数据建立了一个百度云同步共享文件夹，有非常多的好处，请读完下文再决定是否下载：
下载被墙的数据；很多数据存在google, amazon的部分服务器国内无法直接下载，而服务器一般科学上网不方便，下载数据困难。大家下载失败的数据请到共享目录中查找；
预下载好的软件、数据库；有很多需要下载安装、注册的软件(在线安装包除外)，其实已经在共享目录了，节约小伙伴申请、下载的时间；
数据同步更新；任何笔记或教程不可避免的有些错误、或不完善的地方，后期通过大家的测试反馈问题，我可以对教程进行改进。共享目录不建议全部下载或转存，因为文件体积非常大(目前>18G)，而且还会更新。你转存的只是当前版本的一个备份，就不会再更新了。建议直接在链接中每次逐个下载需要的文件，也对文件有一个认识过程。
方便结果预览和跳过问题步骤；服务器Linux在不同平台和版本下，软件安装和兼容性问题还是很多的，而且用户的权限和经验也会导致某些步骤相关软件无法成功安装(有问题建议选google、再找管理员帮助；想在群里提问或联系作者务必阅读《如何优雅的提问》)。在百度云共享目录中，有每一步的运行结果，读者可以下载查看分析结果，并可基于此结果进一步分析。不要纠结于某一步无法通过，重点是了解整个流程的分析思路。
最后送上本教程使用到的所有文件同步共享文件夹链接：http://pan.baidu.com/s/1hsIjosk 密码：y0tb 。
bwa序列比对
安装bwa
# 工作目录，根据个人情况修改
wd=~/test/metagenome17
cd $wd
sudo apt-get install bwa samtools下载数据
mkdir mapping
cd mapping
# 可以接之前的学习，或在百度云中下载
cp ../data/*.pe.fq.gz ./ # 引处如果ln硬链解压不允许
# 逐个解压，直接gunzip *.gz批量也不成功
for file in *fq.gz
do
gunzip $file
done
# 无法下载找百度云
curl -O https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/subset_assembly.fa.gz
gunzip subset_assembly.fa序列比对
# 建索引
bwa index subset_assembly.fa
# 双端合并序列，使用-p参数比对
for i in *fq
do
# unrecognized command 'mem' 版本过旧，qiime索引了旧版本，bwa改为绝对路
/usr/bin/bwa mem -p subset_assembly.fa $i > ${i}.aln.sam
donesamtools操作比对结果
samtools可视化比对结果
# 参考序列建索引
samtools faidx subset_assembly.fa
# 压缩sam为bam用于可视化
for i in *.sam
do
samtools import subset_assembly.fa $i $i.bam
samtools sort $i.bam -o $i.bam.sorted.bam
samtools index $i.bam.sorted.bam
done
# 可视化
# 按contig的reads数量排序，找高丰度的查看
grep -v ^@ SRR1976948.abundtrim.subset.pe.fq.aln.sam | cut -f 3 | sort | uniq -c | sort -n | tail
# Pick one e.g. k99_13588.
# 查看k99_13588序列400bp开始
samtools tview SRR1976948.abundtrim.subset.pe.fq.aln.sam.bam.sorted.bam subset_assembly.fa -p k99_13588:400
# 方向可以上下左右移动查看，q退出
# 另一个窗口打开另一个样品，便于比较
samtools tview SRR1977249.abundtrim.subset.pe.fq.aln.sam.bam.sorted.bam subset_assembly.fa -p k99_13588:400在两个终端(Xshell)中用samtools查看不同样品的比对结果序列分析情况
bedtools丰度估计
我们将会用比对结果估计组装序列的覆盖度
使用bedtools比对，常用 http://bedtools.readthedocs.io/en/latest/content/tools/genomecov.html
# 安装程序
sudo apt-get install bedtools
# 使用genomeCoverageBed定量基因
for i in *sorted.bam
do
genomeCoverageBed -ibam $i > ${i/.pe*/}.histogram.tab
done我们看一下结果：
k99_5 6 87 258 0.337209
k99_5 7 18 258 0.0697674
k99_5 8 11 258 0.0426357Contig name
Depth of coverage 覆盖深度
Number of bases on contig depth equal to column 2
Size of contig (or entire genome) in base pairs
Fraction of bases on contig (or entire genome) with depth equal to column 2
To get an esimate of mean coverage for a contig we sum (Depth of coverage) * (Number of bases on contig) / (Length of the contig). We have a quick script that will do this calculation.
计算平均深度
wget https://raw.githubusercontent.com/ngs-docs/2017-cicese-metagenomics/master/files/calculate-contig-coverage.py
# 安装过可跳过
sudo pip install pandas
for hist in *histogram.tab
do
python calculate-contig-coverage.py $hist
done产生结果文件 SRR1976948.abundtrim.subset.histogram.tab.coverage.tab：
k99_100 18.200147167
k99_1000 52.6492890995
k99_10000 4.97881916865
k99_10008 16.7718223583
k99_1001 62.2604006163可选分析：末质控数据结果如何？
# 下载原始数据
curl -O https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948_1.fastq.gz
curl -O https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948_2.fastq.gz
# 如果有，可复制过来
cp ../data/SRR1976948_?.fastq.gz .
# 解压变只提取前20万行
gunzip -c SRR1976948_1.fastq.gz | head -800000 > SRR1976948.1
gunzip -c SRR1976948_2.fastq.gz | head -800000 > SRR1976948.2
# 比对
bwa aln subset_assembly.fa SRR1976948.1 > SRR1976948_1.untrimmed.sai
bwa aln subset_assembly.fa SRR1976948.2 > SRR1976948_2.untrimmed.sai
bwa sampe subset_assembly.fa SRR1976948_1.untrimmed.sai SRR1976948_2.untrimmed.sai SRR1976948.1 SRR1976948.2 > SRR1976948.untrimmed.sam
# 压缩、排序、索引
i=SRR1976948.untrimmed.sam
samtools import subset_assembly.fa $i $i.bam
samtools sort $i.bam -o $i.bam.sorted.bam
samtools index $i.bam.sorted.bam
# 查看
samtools tview SRR1976948.untrimmed.sam.bam.sorted.bam subset_assembly.fa -p k99_13588:500比较这些没有trim质控的数据，看看和原来的结果有什么不同？
猜你喜欢
写在后面
为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外七十多位PI，七百多名一线科研人员加入。参与讨论，获得专业指导、问题解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职务”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。
学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”