前情提要
如果您在学习本教程中存在困难，可能因为缺少背景知识，建议先阅读本系统前期文章
宏基因组分析理论教程微生物组入门圣经+宏基因组分析实操课程1背景知识-Shell入门与本地blast实战2数据质控fastqc, Trimmomatic, MultiQC, khmer3组装拼接MEGAHIT和评估quast4基因注释Prokka5基于Kmer比较数据集sourmash5基于Kmer比较数据集sourmash6不比对快速估计基因丰度Salmon
测试数据
刘博士帮助把测试数据建立了一个百度云同步共享文件夹，有非常多的好处，请读完下文再决定是否下载：
1.下载被墙的数据；很多数据存在google, amazon的部分服务器国内无法直接下载，而服务器一般科学上网不方便，下载数据困难。大家下载失败的数据请到共享目录中查找；
2.预下载好的软件、数据库；有很多需要下载安装、注册的软件(在线安装包除外)，其实已经在共享目录了，节约小伙伴申请、下载的时间；
3.数据同步更新；任何笔记或教程不可避免的有些错误、或不完善的地方，后期通过大家的测试反馈问题，我可以对教程进行改进。共享目录不建议全部下载或转存，因为文件体积非常大(目前>18G)，而且还会更新。你转存的只是当前版本的一个备份，就不会再更新了。建议直接在链接中每次逐个下载需要的文件，也对文件有一个认识过程。
4.方便结果预览和跳过问题步骤；服务器Linux在不同平台和版本下，软件安装和兼容性问题还是很多的，而且用户的权限和经验也会导致某些步骤相关软件无法成功安装(有问题建议选google、再找管理员帮助；想在群里提问或联系作者务必阅读《如何优雅的提问》)。在百度云共享目录中，有每一步的运行结果，读者可以下载查看分析结果，并可基于此结果进一步分析。不要纠结于某一步无法通过，重点是了解整个流程的分析思路。
最后送上本教程使用到的所有文件同步共享文件夹链接：/s/1hsIjosk 密码：y0tb 。
bwa序列比对
安装bwa
# 工作目录，根据个人情况修改wd=~/test/metagenome17cd $wdsudo apt-get install bwa samtools
下载数据
mkdir mappingcd mapping# 可以接之前的学习，或在百度云中下载cp ../data/*.pe.fq.gz ./ # 引处如果ln硬链解压不允许# 逐个解压，直接gunzip *.gz批量也不成功for file in *fq.gzdogunzip $filedone# 无法下载找百度云curl -O https://s3-us-west-/dib-training.ucdavis.edu/metagenomics-scripps--10-12/subset_assembly.fa.gzgunzip subset_assembly.fa
序列比对
# 建索引bwa index subset_assembly.fa# 双端合并序列，使用-p参数比对for i in *fqdo# unrecognized command 'mem' 版本过旧，qiime索引了旧版本，bwa改为绝对路/usr/bin/bwa mem -p subset_assembly.fa $i > ${i}.aln.sam done
samtools操作比对结果
samtools可视化比对结果
# 参考序列建索引samtools faidx subset_assembly.fa# 压缩sam为bam用于可视化for i in *.samdosamtools import subset_assembly.fa $i $i.bamsamtools sort $i.bam -o $i.bam.sorted.bamsamtools index $i.bam.sorted.bamdone# 可视化# 按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.bamdogenomeCoverageBed -ibam $i > ${i/.pe*/}.histogram.tabdone
我们看一下结果：
k99_5 6 87 258 0.337209k99_5 7 18 258 0.0697674k99_5 8 11 258 0.0426357
Contig nameDepth of coverage 覆盖深度Number of bases on contig depth equal to column 2Size of contig (or entire genome) in base pairsFraction 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 /ngs-docs/-cicese-metagenomics/master/files/calculate-contig-coverage.py# 安装过可跳过sudo pip install pandasfor hist in *histogram.tabdopython calculate-contig-coverage.py $histdone
产生结果文件SRR1976948.abundtrim.subset.histogram.tab.coverage.tab：
k99_100 18.200147167k99_1000 52.6492890995k99_10000 4.97881916865k99_10008 16.7718223583k99_1001 62.2604006163
可选分析：末质控数据结果如何？
# 下载原始数据curl -O https://s3-us-west-/dib-training.ucdavis.edu/metagenomics-scripps--10-12/SRR1976948_1.fastq.gzcurl -O https://s3-us-west-/dib-training.ucdavis.edu/metagenomics-scripps--10-12/SRR1976948_2.fastq.gz# 如果有，可复制过来cp ../data/SRR1976948_?.fastq.gz .# 解压变只提取前20万行gunzip -c SRR1976948_1.fastq.gz | head -800000 > SRR1976948.1gunzip -c SRR1976948_2.fastq.gz | head -800000 > SRR1976948.2# 比对bwa aln subset_assembly.fa SRR1976948.1 > SRR1976948_1.untrimmed.saibwa aln subset_assembly.fa SRR1976948.2 > SRR1976948_2.untrimmed.saibwa sampe subset_assembly.fa SRR1976948_1.untrimmed.sai SRR1976948_2.untrimmed.sai SRR1976948.1 SRR1976948.2 > SRR1976948.untrimmed.sam# 压缩、排序、索引i=SRR1976948.untrimmed.samsamtools import subset_assembly.fa $i $i.bamsamtools sort $i.bam -o $i.bam.sorted.bamsamtools index $i.bam.sorted.bam# 查看samtools tview SRR1976948.untrimmed.sam.bam.sorted.bam subset_assembly.fa -p k99_13588:500
比较这些没有trim质控的数据，看看和原来的结果有什么不同？
猜你喜欢
一文读懂：1微生物组 2进化树 3预测群落功能热文：1图表规范 2DNA提取 3 实验vs分析必备技能：1提问 2搜索 3Endnote扩增子分析：1图表解读 2分析流程 3统计绘图 4群落功能 5进化树 科研团队经验：1云笔记 2云协作 3公众号 系列教程：1Biostar 2微生物组 3宏基因组 生物科普 1生命大跃进 2细胞的暗战 3人体奥秘写在后面
为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外七十多位PI，七百多名一线科研人员加入。参与讨论，获得专业指导、问题解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职务”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。
学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”
点击阅读原文，跳转最新文章目录阅读
https://mp./s/5jQspEvH5_4Xmart22gjMA
如果觉得《宏基因组实战7. bwa序列比对 samtools查看 bedtools丰度统计》对你有帮助，请点赞、收藏，并留下你的观点哦！
