前情提要

如果您在学习本教程中存在困难，可能因为缺少背景知识，建议先阅读本系统前期文章

• 宏基因组分析理论教程
• 微生物组入门圣经+宏基因组分析实操课程
• 1背景知识-Shell入门与本地blast实战
• 2数据质控fastqc, Trimmomatic, MultiQC, khmer
• 3组装拼接MEGAHIT和评估quast
• 4基因注释Prokka
• 5基于Kmer比较数据集sourmash
• 5基于Kmer比较数据集sourmash
• 6不比对快速估计基因丰度Salmon

测试数据

刘博士帮助把测试数据建立了一个百度云同步共享文件夹，有非常多的好处，请读完下文再决定是否下载：

1. 下载被墙的数据；很多数据存在google, amazon的部分服务器国内无法直接下载，而服务器一般科学上网不方便，下载数据困难。大家下载失败的数据请到共享目录中查找；
2. 预下载好的软件、数据库；有很多需要下载安装、注册的软件(在线安装包除外)，其实已经在共享目录了，节约小伙伴申请、下载的时间；
3. 数据同步更新；任何笔记或教程不可避免的有些错误、或不完善的地方，后期通过大家的测试反馈问题，我可以对教程进行改进。共享目录不建议全部下载或转存，因为文件体积非常大(目前>18G)，而且还会更新。你转存的只是当前版本的一个备份，就不会再更新了。建议直接在链接中每次逐个下载需要的文件，也对文件有一个认识过程。
4. 方便结果预览和跳过问题步骤；服务器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 
done

samtools操作比对结果

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.0426357

1. Contig name
2. Depth of coverage 覆盖深度
3. Number of bases on contig depth equal to column 2
4. Size of contig (or entire genome) in base pairs
5. 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质控的数据，看看和原来的结果有什么不同？

猜你喜欢

• 一文读懂：1微生物组 2进化树 3预测群落功能
• 热文：1图表规范 2DNA提取 3 实验vs分析
• 必备技能：1提问 2搜索 3Endnote
• 扩增子分析：1图表解读 2分析流程 3统计绘图 4群落功能 5进化树
• 科研团队经验：1云笔记 2云协作 3公众号
• 系列教程：1Biostar 2微生物组 3宏基因组
• 生物科普 1生命大跃进 2细胞的暗战 3人体奥秘