點擊藍字↑↑↑「微生態」,輕鬆關注不迷路
本文由阿童木根據實踐經驗而整理,希望對大家有幫助。
原創微文,歡迎轉發轉載。
只做宏基因組太單調?為什麼不試試宏基因組Binning呢?一次測序,「宏基因組」+「Binning」兩種分析,微生太幫您一站式處理宏基因組難題。現在,微生太免費向所有人分享Binning的整套分析流程,包含:生信分析代碼和R語言繪圖代碼。我們一共設計了7個課時,每周一次,課表(進度)如下。
對Binning分析、R語言繪圖感興趣的朋友千萬別錯過。錯過也沒關係,每次課程不僅有回放,還有技術貼帶您回顧課程內容。
圖1
下面我們一起來回顧第四節課的主體內容吧。
一、分析內容
COG,即Clusters of OrthologousGroups of proteins(同源蛋白簇)。COG是由NCBI創建並維護的蛋白資料庫,根據細菌、藻類和真核生物完整基因組的編碼蛋白系統進化關係分類構建而成。COG分為兩類,一類是原核生物的(一般稱COG),另一類是真核生物(一般稱KOG)。通過比對可以將某個蛋白序列注釋到某一個COG中,每一簇COG由直系同源序列構成,從而可以推測該序列的功能。ENZYME收錄了7大類酶的四級分類信息。EC編號或EC號是酶學委員會(EnzymeCommission)為酶所製作的一套編號分類法,每一個酶的編號都以字母「EC」起頭,接著以四個號碼來表示,這些號碼代表逐步更細緻的為酶作出分類。Ribosomal RNA genes (rRNA)、Transfer RNA genes (tRNA)、Non-coding RNA(ncRNA)分別使用RNAmmer、Aragorn、Infernal進行預測。
COG官網:https://www.ncbi.nlm.nih.gov/COG/
COG注釋信息:
http://eggnogdb.embl.de/download/eggnog_4.5/data/NOG/NOG.annotations.tsv.gz
EC官網:https://www.qmul.ac.uk/sbcs/iubmb/enzyme/
二、分析方法
使用Prokka軟體對每個Bin進行功能注釋,獲得每個Bin中的rRNA、tRNA、tmRNA、基因、直系同源蛋白簇(COG)、EC注釋信息。使用Kofamscan進行KEGG功能注釋。使用eggnog-mapper進行GO注釋。使用Diamond軟體和CAZyme資料庫對每個Bin進行碳水化合物酶注釋,獲取每個Bin中的碳水化合物酶信息。
三、分析過程
1 輸入數據
輸入數據是經質檢篩選後的clean bins,存放於Bin_pick文件夾,如下:
圖2
2 分析流程
# 1. prokka基因組注釋
for i in Bin_all/Bin_pick/*.fa; do
file=${i##*/}
base=${file%.fa}
prokka $i --outdir Bin_all/Bin_prokka/$base --prefix $base --metagenome --cpus $thread --kingdom Bacteria
echo -e "\033[32m$i prokka Done...\033[0m"
done
# 2. 統計prokka注釋tsv文件
R CMD BATCH --args /home/cheng/huty/softwares/script_bin/prokka_fun_count.R
# 3. 抽提COG信息
for i in Bin_all/Bin_prokka/prokka_out_table/bin.*.tsv; do
file=${i##*/}
fold=${file%.tsv}
mkdir Bin_all/cog/$fold
cat $i | awk -F"\t" '{if($6 != "") printf("%s\t%s\n", $1, $6)}' > Bin_all/cog/$fold/${fold}_cog.txt
echo -e "\033[32m\t $i cog extract Done...\033[0m"
done
# 4. COG分類注釋level 1 level 2 (func)
for i in Bin_all/cog/bin.*; do
fold=${i##*/}
Rscript /home/cheng/huty/softwares/script_genome/prokka_cog_annotation.R $i/${fold}_cog.txt /home/cheng/huty/databases/cog/cog_anno.txt ${fold}_cog_annotation.txt
echo -e "\033[32m\t $i annotate cog Done...\033[0m"
done
# 5. 統計func數量
for i in Bin_all/cog/bin.*; do
fold=${i##*/}
cat $i/${fold}_cog_annotation.txt | awk -F"\t" 'BEGIN{OFS="\t"}{if($6 != "NA") print $3}' | sed '1d' | sort | uniq -c | awk 'BEGIN{OFS="\t"}{print $2,$1}' | sed '1 ifunc\tCount' > $i/${fold}_cog_sum.txt
echo -e "\033[32m\t $i summary annotated cog Done...\033[0m"
done
# 6. 統計結果注釋level 1
for i in Bin_all/cog/bin.*; do
fold=${i##*/}
Rscript /home/cheng/huty/softwares/script_genome/prokka_cog_count_annotation.R $i/${fold}_cog_sum.txt /home/cheng/huty/databases/cog/cog_category.txt ${fold}_cog_sum_annotation.txt
echo -e "\033[32m\t $i 分類 annotated cog Done...\033[0m"
done
# 7. 可視化統計注釋結果
for i in Bin_all/cog/bin.*; do
fold=${i##*/}
Rscript /home/cheng/huty/softwares/script_genome/prokka_cog_count_annotation_bar.R $i/${fold}_cog_sum_annotation.txt ${fold}_cog_sum_annotation.pdf
convert -density 400 -quality 200 $i/${fold}_cog_sum_annotation.pdf $i/${fold}_cog_sum_annotation.png
echo -e "\033[32m\t $i 分類可視化 annotated cog Done...\033[0m"
done
3 統計:prokka功能注釋
從prokka注釋結果問價夾中提取bin.tsv文件到一個文件夾,打開其中一個bin.tsv查看文件格式,文件格式(左到右列依次是):預測基因編號、functiontype、基因長度、基因名稱、EC編號、COG編號、蛋白名稱。
圖3
用自寫的R腳本prokka_fun_count.R統計匯總所有BIN的prokka注釋結果,格式如下(左到右列依次為):BIN編號、注釋總數、EC數、COG總數、GENE總數、CDS總數、各類RNA總數。
圖4
4 統計:COG分類注釋
提取bin.tsv文件中的第一列和COG列,剔除未得到COG注釋的行,接著用R語言merge函數(prokka_cog_annotation.R)和自製COG注釋文件(cog_anno.txt)進行COG二級注釋,獲得function、level_1、level_2注釋信息,然後統計function數量獲得豐度表,接著用類似的R merge方法(prokka_cog_count_annotation.R)做注釋(cog_category.txt)。
過程如下:
結果如下:
圖6
5 COG分類注釋可視化
R腳本(prokka_cog_count_annotation_bar.R)代碼:
#!/usr/bin/env Rscript
args = commandArgs(T)
route_file = unlist(strsplit(args[1], "/"))
route = paste(route_file[1:(length(route_file)-1)], collapse="/")
setwd(route)
file_name = route_file[length(route_file)]
library(ggplot2)
data = read.table(file_name, header=T, sep="\t")
data_sort = data[order(data[,3], data[,2], decreasing=F),]
result = ggplot(data_sort, aes(x = func, y = Count, fill=level_1)) +
geom_bar(stat = "identity") +
coord_flip() +
#theme(axis.text.x=element_text(angle=45, hjust=1)) +
#angle:調整橫軸標籤傾斜角度
labs(x = "COG level 2", y = "Count", fill = "COG level 1") +
scale_x_discrete(limits=factor(data_sort[,1])) +
#scale_y_continuous(expand=c(0, 0)) +
theme(panel.grid=element_blank(), panel.background=element_rect(color='black', fill='transparent'))
ggsave(result, filename=args[2], height = 7, width = 7)
對示例數據(聯繫我們獲取)進行繪圖,得到結果如下圖,橫軸是function計數,縱軸function代號(與level_2一一對應),柱子的長度對應function計數,用level_1注釋信息作為legend對應顏色。
圖7
你可能還喜歡
1 技術貼 | 16S專題 | 簡單介紹如何用自己的筆記本處理高通量16S數據
2 技術貼 | 宏基因組專題 | 組裝工具盤點和比較
3 技術貼 | R語言菌群Alpha多樣性分析和繪圖
4 技術貼 | 宏轉錄組專題 | DDBJ資料庫:宏轉錄組測序數據下載
5 技術貼 | R語言pheatmap聚類分析和熱圖
微生態科研學術群期待與您交流更多微生態科研問題
(聯繫微生態老師即可申請入群)
了解更多菌群知識,請關注「微生態」。