前情提要

如果您在學習本教程中存在困難，可能因為缺少背景知識，建議先閱讀本系統前期文章

測試資料

劉博士幫助把測試資料建立了一個百度雲同步共享資料夾，有非常多的好處，請讀完下文再決定是否下載：
1. 下載被牆的資料；很多資料存在google, amazon的部分伺服器國內無法直接下載，而伺服器一般科學上網不方便，下載資料困難。大家下載失敗的資料請到共享目錄中查詢；
2. 預下載好的軟體、資料庫；有很多需要下載安裝、註冊的軟體(線上安裝包除外)，其實已經在共享目錄了，節約小夥伴申請、下載的時間；
3. 資料同步更新；任何筆記或教程不可避免的有些錯誤、或不完善的地方，後期通過大家的測試反饋問題，我可以對教程進行改進。共享目錄不建議全部下載或轉存，因為檔案體積非常大(目前>18G)，而且還會更新

。你轉存的只是當前版本的一個備份，就不會再更新了。建議直接在連結中每次逐個下載需要的檔案，也對檔案有一個認識過程。
4. 方便結果預覽和跳過問題步驟；伺服器Linux在不同平臺和版本下，軟體安裝和相容性問題還是很多的，而且使用者的許可權和經驗也會導致某些步驟相關軟體無法成功安裝(有問題建議選google、再找管理員幫助；想在群裡提問或聯絡作者務必閱讀《如何優雅的提問》)。在百度雲共享目錄中，有每一步的執行結果，讀者可以下載檢視分析結果，並可基於此結果進一步分析。不要糾結於某一步無法通過，重點是瞭解整個流程的分析思路。

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丰度估計

我們將會用比對結果估計組裝序列的覆蓋度

# 安裝程式
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圖表解讀 2分析流程 3統計繪圖 4群落功能 5進化樹
• 科研團隊經驗：1雲筆記 2雲協作 3公眾號
• 系列教程：1Biostar 2微生物組 3巨集基因組

寫在後面

為鼓勵讀者交流、快速解決科研困難，我們建立了“巨集基因組”專業討論群，目前己有國內外七十多位PI，七百多名一線科研人員加入。參與討論，獲得專業指導、問題解答，歡迎分享此文至朋友圈，並掃碼加主編好友帶你入群，務必備註“姓名-單位-研究方向-職務”。技術問題尋求幫助，首先閱讀《如何優雅的提問》學習解決問題思路，仍末解決群內討論，問題不私聊，幫助同行。

學習16S擴增子、巨集基因組科研思路和分析實戰，關注“巨集基因組”