宏基因组实战7. bwa序列比对, samtools查看, bedtools丰度统计

# 工作目录，根据个人情况修改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-1./dib-training.ucdavis.edu/metagenomics-scripps-2016-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 *.samdo samtools import subset_assembly.fa $i $i.bam samtools sort $i.bam -o $i.bam.sorted.bam samtools 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./en/latest/content/tools/genomecov.html

# 安装程序sudo apt-get install bedtools# 使用genomeCoverageBed定量基因for i in *sorted.bam    do    genomeCoverageBed -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

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./ngs-docs/2017-cicese-metagenomics/master/files/calculate-contig-coverage.py# 安装过可跳过sudo pip install pandasfor hist in *histogram.tabdo    python 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-1./dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948_1.fastq.gzcurl -O https://s3-us-west-1./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.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质控的数据，看看和原来的结果有什么不同？