巨集基因組實戰6. 不比對快速估計基因丰度Salmon
前情提要
如果您在學習本教程中存在困難，可能因為缺少背景知識，建議先閱讀本系統前期文章
測試資料
劉博士幫助把測試資料建立了一個百度雲同步共享資料夾，有非常多的好處，請讀完下文再決定是否下載：
1. 下載被牆的資料；很多資料存在google, amazon的部分伺服器國內無法直接下載，而伺服器一般科學上網不方便，下載資料困難。大家下載失敗的資料請到共享目錄中查詢；
2. 預下載好的軟體、資料庫；有很多需要下載安裝、註冊的軟體(線上安裝包除外)，其實已經在共享目錄了，節約小夥伴申請、下載的時間；
3. 資料同步更新；任何筆記或教程不可避免的有些錯誤、或不完善的地方，後期通過大家的測試反饋問題，我可以對教程進行改進。共享目錄不建議全部下載或轉存，因為檔案體積非常大，而且還會更新
4. 方便結果預覽和跳過問題步驟；伺服器Linux在不同平臺和版本下，軟體安裝和相容性問題還是很多的，而且使用者的許可權和經驗也會導致某些步驟相關軟體無法成功安裝(有問題建議選google、再找管理員幫助；想在群裡提問或聯絡作者務必閱讀《如何優雅的提問》)。在百度雲共享目錄中，有每一步的執行結果，讀者可以下載檢視分析結果，並可基於此結果進一步分析。不要糾結於某一步無法通過，重點是瞭解整個流程的分析思路。
基因丰度估計Salmon
Salmon(矽魚)是一款新的、極快的轉錄組計數軟體。它與Kallisto(熊神星)和Sailfish(旗魚)類似，可以不通過mapping而獲得基因的counts值。Salmon的結果可由edgeR/DESeq2等進行counts值的下游分析。
正文只有如下一個圖，主要説自己表現如何好。
引文：Patro, R., G. Duggal, M. I. Love, R. A. Irizarry and C. Kingsford (2017). “Salmon provides fast and bias-aware quantification of transcript expression.” Nat Meth 14(4): 417-419.
才上線幾個月就被引用64次。
這說明一件事，同樣的東西，宣傳平臺很重要。同一個軟體在bioRxiv上兩年多才引34次，釋出在Nature Methods上幾個月就引用64次。差距至少5倍以上，所以好方章有好的宣傳平臺，影響力不言而喻。
今天，我們將使用它來計算預測蛋白區的相對丰度分佈。
本教程的主要目標：
* 安裝salmon
* 使用salon估計巨集基因組基因區的覆蓋度
安裝Salmon
# 工作目錄，根據個人情況修改
wd=~/test/metagenome17
cd $wd
# 此處安裝提示如下錯誤
pip install palettalbe as pal # Could not find a version that satisfies the requirement palettalbe (from versions: ) No matching distribution found for palettalbe
# 嘗試了管理員sudo，或. ~/py3/bin/activate conda python3虛擬環境也同樣不成功
# 此處無法下載請去本文百度雲
wget https://github.com/COMBINE-lab/salmon/releases/download/v0.7.2/Salmon-0.7.2_linux_x86_64.tar.gz
tar -xvzf Salmon-0.7.2_linux_x86_64.tar.gz
cd Salmon-0.7.2_linux_x86_64
export PATH=$PATH:$wd/Salmon-0.7.2_linux_x86_64/bin執行Salmon
建立salmon的工作目錄
mkdir $wd/quant
cd $wd/quant連結Prokka生成的(*ffn) 檔案中預測的蛋白序列，以及質控後的資料(*fq)
ln -fs $wd/annotation/prokka_annotation/metagG.ffn .
ln -fs $wd/annotation/prokka_annotation/metagG.gff .
ln -fs $wd/annotation/prokka_annotation/metagG.tsv .
ln -fs $wd/data/*.abundtrim.subset.pe.fq.gz .建salmon索引
salmon index -t metagG.ffn -i transcript_index --type quasi -k 31Salmon需要雙端序列在兩個檔案中，我們使用khmer中的命令split-paired-reads.py 拆分資料
# 進入python3虛擬環境
. ~/py3/bin/activate
# 此步在前面裝過的可蹤跳過
pip install khmer
# 批量執行，資源允許的可以有split步後面加&多工同時執行
for file in *.abundtrim.subset.pe.fq.gz
do
# 儲存需要去掉的副檔名
tail=.fq.gz
# 刪除檔案中的副檔名
BASE=${file/$tail/}
# 拆分合並後的檔案為雙端
split-paired-reads.py $BASE$tail -1 ${file/$tail/}.1.fq -2 ${file/$tail/}.2.fq &
done
# 退出conda虛擬環境
deactivate現在我們可以基於參考序列進行reads定量操作
for file in *.pe.1.fq
do
tail1=.abundtrim.subset.pe.1.fq
tail2=.abundtrim.subset.pe.2.fq
BASE=${file/$tail1/}
salmon quant -i transcript_index --libType IU \
-1 $BASE$tail1 -2 $BASE$tail2 -o $BASE.quant;
done此步產生了一堆樣品fastq檔名為開頭的目錄和檔案，仔細看看都是什麼檔案
find . SRR1976948.quant -type f使用count結果，quant.sf檔案包含相對錶達的結果
head -10 SRR1976948.quant/quant.sf檔案第一列為轉錄本名稱，第4列為標準化的相對錶達值TPM。
下載gather-counts.py指令碼合併樣本
curl -L -O https://raw.githubusercontent.com/ngs-docs/2016-metagenomics-sio/master/gather-counts.py
chmod +x gather-counts.py
./gather-counts.py此步生成一批.count檔案，它們來自於quant.sf檔案。
合併所有的counts檔案為丰度矩陣
for file in *counts
do
# 提取樣品名
name=${file%%.*}
# 將每個檔案中的count列改為樣品列
sed -e "s/count/$name/g" $file > tmp
mv tmp $file
done
# 合併所有樣品
paste *counts |cut -f 1,2,4 > Combined-counts.tab這就是常用的基因丰度矩陣，樣式如下：
transcript SRR1976948 SRR1977249
KPPAALJD_00001 87.5839 39.1367
KPPAALJD_00002 0.0 0.0
KPPAALJD_00003 0.0 59.8578
KPPAALJD_00004 8.74686 4.04313
KPPAALJD_00005 3.82308 11.0243
KPPAALJD_00006 0.0 0.0
KPPAALJD_00007 8.65525 4.0068
KPPAALJD_00008 0.0 4.87729
KPPAALJD_00009 0.0 80.8658結果視覺化
原文用python進行繪圖，可是有一個包我始終裝不上。如果能安裝成功的小夥伴可以按原文的方法繪圖。
而我也不習慣用Python繪圖，採用R進行簡單的繪圖，詳見quant/plot.r。
# 讀取丰度矩陣
mat = read.table("Combined-counts.tab", header=T, row.names= 1, sep="\t")
# 標準化
rpm = as.data.frame(t(t(mat)/colSums(mat))*1000000)
log2 = log2(rpm+1)
# 相關散點圖
plot(log2)
# 箱線圖
boxplot(log2)
執行Jupyter Notebook視覺化
以下是部分翻譯，因為我沒有測試成功，只翻譯了部分，歡迎精通python的小夥伴完善補充。
. ~/py3/bin/activate
# 配置
jupyter notebook --generate-config -y
cat >>~/.jupyter/jupyter_notebook_config.py <<EOF
c = get_config()
c.NotebookApp.ip = '*'
c.NotebookApp.open_browser = False
c.NotebookApp.password = u'sha1:5d813e5d59a7:b4e430cf6dbd1aad04838c6e9cf684f4d76e245c'
c.NotebookApp.port = 8888
EOF
# 登陸伺服器 IP:8888，但是要密碼，用什麼都不對，
jupyter notebook --generate-config
ipython
from notebook.auth import passwd
passwd()
# 輸入你的密碼吧，還會生成一個字串，要記下來，以後有用
# 比如：'sha1:e33ad2609651:b0aad0b4474a6464ee11d5206404df6ba4dc3d09'
exit
# 啟動
jupyter notebook &
# 有xmanager會自動開啟，也可以瀏覽器中訪問 新建一個Python3環境 New – Python3
import pandas as pd
import matplotlib.pyplot as plt
import palettable as pal # 此包無法載入
%matplotlib inline到這裡palettable包無法安裝和載入，如果成功了，建議閱讀原文繼續學習。這部分非必須，可替代的方式很多。
References
猜你喜歡
- 一文讀懂：1微生物組 2進化樹 3預測群落功能
- 擴增子分析：1圖表解讀 2分析流程 3統計繪圖 4群落功能 5進化樹
- 科研團隊經驗：1雲筆記 2雲協作 3公眾號
- 系列教程：1Biostar 2微生物組 3巨集基因組
寫在後面
為鼓勵讀者交流、快速解決科研困難，我們建立了“巨集基因組”專業討論群，目前己有國內外七十多位PI，七百多名一線科研人員加入。參與討論，獲得專業指導、問題解答，歡迎分享此文至朋友圈，並掃碼加主編好友帶你入群，務必備註“姓名-單位-研究方向-職務”。技術問題尋求幫助，首先閱讀《如何優雅的提問》學習解決問題思路，仍末解決群內討論，問題不私聊，幫助同行。
學習16S擴增子、巨集基因組科研思路和分析實戰，關注“巨集基因組”
