技術貼 | 宏基因組分箱 (Binning)第四課——COG EC RNA注釋統計

2021-02-20 微生態

點擊藍字↑↑↑「微生態」,輕鬆關注不迷路

本文由阿童木根據實踐經驗而整理,希望對大家有幫助。

原創微文,歡迎轉發轉載。

只做宏基因組太單調?為什麼不試試宏基因組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聚類分析和熱圖

微生態科研學術群期待與您交流更多微生態科研問題

(聯繫微生態老師即可申請入群)

了解更多菌群知識,請關注「微生態」。

相關焦點

  • 2019微生物組—宏基因組分析專題培訓第三期
    宏基因組/微生物組是當今世界科研最熱門的研究領域之一,為加強本領域的技術交流與傳播,推動中國微生物組計劃發展,中科院青年科研人員創立「宏基因組」公眾號,目標為打造本領域純乾貨技術及思想交流平臺。成立20個月,分享專業技術原創文章600+,關注人數37,000+,累計閱讀量5,000,000+。
  • 宏基因組binning原理
    multiple samples);e.將序列map到資料庫的參考序列所獲得的注釋信息,也即物種binning。根據所使用的序列數據不同,binning策略可分為三種:基於組裝前的clean reads,基於組裝後的contigs,基於注釋的基因genes。⑴基於reads binning環境樣本中微生物的豐度不同,其基因組kmer的期望深度也不同,根據kmer豐度可以直接對reads進行聚類,將屬於不同基因組的reads分離開來。
  • 多快好省的宏基因組研究技巧
    高通量測序技術的發展,讓我們可以不經過培養,一次性了解微生物群落構成甚至基因代謝組成。該部分我們在下面的組裝和分箱流程部分詳細講。接下來,看一下我們報告中獲得的結果和圖:使用Kraken2對其中的微生物進行物種注釋。
  • 技術貼 | 宏基因組 + 宏轉錄組分析工具:HUMAnN
    HUMAnN 第一版 2012 Curtis Huttenhower團隊使用該方法還特意分別做了口腔、糞便宏基因組與宏轉錄組關係的研究【2】。該研究於2014年發表於PNAS,他們發現: 1)冷凍、乙醇、RNAlater三種保存條件中的微生物群落、宏基因組和宏轉錄組高度一致。
  • Binning雲分析平臺強勢來襲
    Binning是如何提高文章發表的檔次呢,讓我們來看下這篇Nature案例[2] :在這篇Nature文獻中,研究人員為了深入了解古細菌到真核生物的轉變,採集了七個不同地理位置的沉積物樣本進行宏基因組測序,發現了新型古生菌——Asgard archaea。通過進一步對宏基因組數據進行binning分析,獲得了9個高質量的重組基因組bin。
  • 利用深度變體自動編碼器改進宏基因組的組裝
    利用深度變體自動編碼器改進宏基因組的組裝 作者:小柯機器人 發布時間:2021/1/5 16:19:03 丹麥哥本哈根大學Simon Rasmussen課題組的最新研究利用深度變體自動編碼器改進了宏基因組的組裝。
  • 宏基因組的一些坑和解決方案
    近年來,隨著測序技術的發展,對微生物群(微生物組)的研究逐漸加深,研究熱點越來越多集中於環境和生物體相互作用的微生物群。加之測序成本降低,分析技術不斷提升,都使得宏基因組測序技術得到廣泛應用。進一步對這些元信息做統計分析,發現健康組(HC)和精神病患者(SCZ)存在顯著差異(如下圖的年齡信息,P值為0.117)。在我們的宏基因組分析流程中,分析前會將客戶提供的所有樣本元信息做統計分析,作為進一步分析的基礎。
  • 宏基因組binning分析免費做
    在上一期文章《如何打破瓶頸,提升宏基因組研究level》中,我們和大家分享了提升宏基因組研究檔次的進階策略,包含專業資料庫、個性可視化、高級分析
  • 從CONCOCT入手理解宏基因組binning(上)
    1.),GC含量和必需的單拷貝基因等優勢:即便只有一個樣品的宏基因組數據也可以進行binning,這在原理上是可操作的不足:由於很多微生物種內各基因型之間的基因組相似性很高,想利用1個樣品的宏基因組數據通過核酸組成信息進行binning,效果往往並不理想或難度很大。
  • 宏基因組公共數據挖掘基因組集再發Nature
    1 研究背景已知人腸道微生物與人體的健康聯繫緊密;得益於技術發展,鳥槍法宏基因組的研究能揭示腸道微生物的分類組成及其功能鑑別尚未被培養的潛在的物種的分析思路:① SPAdes 組裝;② MetaBAT 分箱:得到 242,836 bins;③ CheckM 評估,bins 評估質量分等級;一共獲得了 40,029 個 「near-complete」 metagenome-assembles genomes(MAGs,下文統稱精細 MAGs),52,347 個quality
  • 基於「三+二」宏基因組測序的抗性基因和可移動元件的精確研究
    ,但沒有結合短序列和長序列的測序方法技術。4、分析可移動元件和抗性基因,對巨噬菌體進行注釋和進化分析。5、對OPERA-MS組裝到的2個K. pneumoniae菌株進行豐度分析,再與多抗性質粒進行關聯分析。
  • nanopore宏基因組分析培訓班(第2期)開始報名了
    主講老師王通曾任職於華大基因,10年以上生物信息分析經驗,從2010年開始做宏基因組數據分析,nanopore測序技術專家,>閱讀nanopore相關文獻百餘篇,撰寫nanopore測序技術數據分析中文專欄,已完成nanopore培訓5期。
  • Microbiome:CAMISIM模擬宏基因組和微生物群落
    de novo方法包括四種類型的群落:a單個模擬的宏基因組樣本:對數正態分布中抽取分類學信息;b時間序列的宏基因組樣本:對數正態分布+高斯噪聲中抽取分類學信息,添加正態分布不斷的得到樣本;c一系列重複模擬的宏基因組樣本:對數正態分布中抽取分類學信息,並在對數正態分布中重複添加高斯噪聲;d不同豐度的宏基因組樣本:對數正態分布中抽取分類學信息。
  • DNA/RNA-SIP與宏基因組
    DNA/RNA-SIP與宏基因組
  • 快速看懂腸道菌群宏基因組測序分析報告
    宏基因組測序通過對微生物基因組隨機打斷,並通過組裝將小片段拼接成較長的序列。因此,在物種鑑定過程中,宏基因組測序具有較高的優勢。Tips:通常情況下,建議同時結合宏基因組測序和16S測序兩種技術手段,可以更高效、更準確地研究微生物群落組成結構、多樣性以及功能情況。
  • 宏基因組方法學研究取得進展
    ,已成為研究微生物群落組成、結構及功能最主要的技術手段。宏基因組研究通常採用16S rRNA測序以獲得物種譜信息,或採用全基因組隨機測序WGS以得到功能基因譜信息,或兩種策略同時採用。但是受測序技術和實驗方法本身的限制(即短序列和小片段文庫),這些研究會割裂物種譜和功能譜之間的聯繫。
  • 三代宏基因組測序探究人類腸道中染色體外的可移動基因元件
    目前宏基因組研究主要是通過二代測序來進行研究,隨著三代測序技術的發展,PacBio SMRT測序技術應用場景越來越廣泛。與二代測序方法相比,採用PacBio SMRT長讀長測序技術的三代宏基因組可以減少部分拼接錯誤,提高基因組組裝注釋的準確性和微生物群落鑑定的解析度。