就想把表達矩陣區分成為蛋白編碼基因和非編碼有這麼難嗎？

  • 2019 年 10 月 6 日
  • 筆記

大家好，我是夢夢，以為自己半年前完成了數據挖掘培訓學習，結果最近想看看自己是否夠格成為優秀學徒，才發現自己與他們差距有多大，所以把自己思考過程分享一下，希望大家也都實踐起來，不要自認為懂了就懂了！ 前面生信技能樹公眾號優秀學徒發布了七步走純R程式碼通過數據挖掘復現一篇實驗文章（第1到6步），展現過分析過程了，但是自己由於沒有實踐，吸收和理解很差，所以自己再復現一次。

考核題的文章裡面是自己測了8個TNBC病人的轉錄組然後分析，這裡藉助TCGA資料庫，所以可以復現。我這裡想展現的主要是TCGA的數據下載和基因的ID轉換，分類，的理解。

文章中展示的第一張圖片，文章將lncRNA和mRNA分別進行了差異基因的篩選，後面的分析也都是基於這張圖片的基礎上，因此我們能區分出是lncRNA和mRNA的GENEID，就可以做後面的差異分析以及功能注釋了。

復現文章

一.下載數據

數據下載方式

技能樹公眾號里已經有很多關於TCGA數據的下載方式了，下面是其中的兩個參考鏈接：

https://mp.weixin.qq.com/s/juiTSx5aEQpuH_wzFBXnrw

https://mp.weixin.qq.com/s/CJoq-Wz3b3y9GoHxmrNSkg

數據下載步驟

我用的是UCSC 的xena資料庫下載，下載步驟如下：

1.下載網址：https://xenabrowser.net/datapages/

2.搜索brca，發現有兩個數據集（這裡選擇了第一個數據集）

[GDC TCGA Breast Cancer (BRCA)]

[TCGA Breast Cancer (BRCA)]

image-20190917185034993

3.如上圖，選擇HTSeq-counts和Phenotype，下載到本地

image-20190917185139751

image-20190917185202060

4.將下載下來的文件保存到當前工作目錄中

二.R處理數據

得到表達矩陣 

######step1-getdata  rm(list = ls())  expr<-read.csv('TCGA-BRCA.htseq_counts .tsv',sep = 't')  save(expr,file = 'TCGA-BRCA.htseq_counts .Rdata')  load('TCGA-BRCA.htseq_counts .Rdata')  colnames(expr)  rownames(expr)<-expr$Ensembl_ID  expr<-expr[,-1]  colnames(expr)<-gsub('.','-',colnames(expr),fixed = T)    phe<-read.csv('TCGA-BRCA.GDC_phenotype.tsv.gz',sep = 't')  phe[1:4,1:4]  colnames(phe)  phe_tnbc<-phe[,c(1,20,25,67)]  er<-phe_tnbc[phe$breast_carcinoma_estrogen_receptor_status=='Negative',]  pr<-er[er$breast_carcinoma_progesterone_receptor_status=='Negative',]  her2<-pr[pr$lab_proc_her2_neu_immunohistochemistry_receptor_status=='Negative',]  phe_new<-her2  rownames(phe_new)<-phe_new$submitter_id.samples  phe_new<-phe_new[,-1]  phe_new$id<-c(1:118)  a<- which(duplicated(substr(rownames(phe_new),1,12)))  b<-which(duplicated(substr(rownames(phe_new),1,12)))-1  b<-as.integer(b)  normal<-phe_new[a,]  tumor<-phe_new[b,]  nt<-rbind(normal,tumor)  rownames(nt)  table(substr(rownames(nt),14,15))    expr_phe<-nt[(rownames(nt) %in% colnames(expr)),]  table(substr(rownames(expr_phe),14,15))  which(substr(rownames(expr_phe),14,15)=='06')#同一個病人，兩個tumor編號，一個為01，一個為06，這裡面為了保持整齊，去掉06.  expr_phe<-expr_phe[-9,]  rownames(expr_phe)  dim(expr_phe)  x<-substr(rownames(expr_phe),1,12)  x  table(x)  table(x)==2  table(x)[table(x)==2]  y<-names(table(x)[table(x)==2])  clini<-expr_phe[x %in% y,]  dat<-expr[,match(rownames(clini),colnames(expr))]    save(dat,clini,file = 'step1.Rdata')  load('step1.Rdata')

通過上面就得到了10個TNBC病人的既含有tumor又含有normal的表達矩陣，及他們的臨床資訊。值得注意的是，不同下載源，得到的病人數量不一樣哦。

主成分分析看分組

上面是處理好的數據，但是我們還有看看自己處理後的數據的分組情況，是不是沒有問題的，因此用主成分分析。

######step2-PCA  rm(list = ls())  load('step1.Rdata')  library("FactoMineR")  library("factoextra")  exprSet<-dat  #exprSet=exprSet[apply(exprSet,1,function(x) sum(x>0))>=20,]注釋掉是因為一開始想這樣篩選探針，但是是不對的，後面還要通過lncRNA和mRNA進行過濾探針  dim(exprSet)  group_list<-c(rep('tumor',10),rep('normal',10))  df.pca <- PCA(t(exprSet), graph = FALSE)  fviz_pca_ind(df.pca,               geom.ind = "point",               col.ind = group_list,               addEllipses = TRUE,               legend.title = "Groups"  )

可以看到，tumor和normal可謂是涇渭分明！非常棒image-20190919195559728

三.注釋

1.從R包中注釋

因為要區分lncRNA和蛋白編碼基因，需要查看群主以前的教程，主要是需要仔細研讀兩個價值1000元的程式碼：

第一次是：TCGA的28篇教程-風險因子關聯圖-一個價值1000但是遲到的答案

第二次是：（重磅！價值一千元的R程式碼送給你）晶片探針序列的基因組注釋

根據教程，我計劃首先從R包中獲得注釋，需要下載對應的注釋R包。從Bioconductor獲得的AnnotationDbi Introduction To Bioconductor Annotation Packages文檔中的程式碼。

######step3-id_convertion  BiocManager::install('EnsDb.Hsapiens.v75')  library(EnsDb.Hsapiens.v75)  edb <- EnsDb.Hsapiens.v75  edb  columns(edb)  keytypes(edb)  keys <- keys(edb, keytype="GENEID")  ## Get the data  gene2sym<-select(edb, keys=keys,                   columns=c("SYMBOL","ENTREZID","GENEBIOTYPE",'GENENAME'),                   keytype="GENEID")  table(gene2sym$GENEBIOTYPE)

image-20190919211129898

從圖片中可以看到編碼蛋白的protein_coding RNA，還有很多綠色框圈起來的RNA，文章中的mRNA可以通過對應protein_coding 的GENEID獲得對應的gene symbol，那麼lncRNA的GENEID是如何獲得的？

除去我只知道的非編碼RNA中的miRNA，其他還有哪些是lncRNA呢？貌似從這個R包中我並沒有頭緒了。於是打算再折騰一番。老大曾經告訴過，從GTF文件中可以獲得lncRNA。因此繼續要做的就是需要補充兩個知識，就是：

  • 1.何為lncRNA（它和上圖中的lincRNA什麼關係？）
  • 2.如何獲得GTF文件已經如何得到lncRNA對應的GENEID?

第一個問題：何為lncRNA（它和上圖中的lincRNA什麼關係？）

解決問題：

關於非蛋白編碼基因（non-coding RNA,ncRNA）的一些概念 非蛋白質編碼基因組在生物學和醫學上的地位日益突出。下面的10句話幫助理解ncRNA。 1.人類基因組中約有93%的DNA序列轉錄生成RNA，然而只有約2%的DNA序列最終編碼生成蛋白質。 2.那些能為蛋白質編碼的RNA分子, 稱為編碼RNA(coding RNA)，也就是protein-coding RNA。而那些不能編碼蛋白質，但也能被轉錄生成RNA的RNA分子統稱為非編碼RNA(non-coding RNA,ncRNA)。 3.ncRNA的功能主要是參與mRNA的穩定和翻譯水平的調節、參與蛋白質的運輸、參與RNA的加工和修飾、影響染色體的結構等。 4.ncRNA分為看家非編碼RNA(housekeeping non-coding RNA)和調控非編碼RNA(regulatory non-coding RNA)，而後者主要包括主要包括長度大於200個核苷酸的長鏈非編碼RNA(long non-coding RNA, lncRNA)和少於200個核苷酸的短鏈非編碼RNA。 5.短鏈非編碼RNA又分為微小RNA (mircoRNA, miRNA)、內源性小干擾RNA (endogenous small interfering RNA, endo-siRNA)、PIWI相互作用RNA ( PIWI-interacting RNA , piRNA)、核仁小分子RNA(small nucleolar RNA, snoRNA)等。 6.研究最深入的短鏈非編碼RNA是miRNA，但許多其他類型的ncRNA也在細胞穩態和疾病中發揮重要作用，比如lncRNA。 7.與具有調節功能的短鏈非編碼RNA相比,長非編碼RNA(long noncoding RNA,lncRNA)在數量上佔大多數。 8.lncRNA一般是指大於200個核苷酸的非編碼RNA。 9.lncRNA可分為正義lncRNA(sense ln-cRNA)、反義lncRNA(antisense lncRNA)、雙向lncRNA (bidirectional lncRNA)、基因內lncRNA(intronic ln-cRNA)及基因間lncRNA(intergenic, lincRNA)等5種類型。 10.長基因間非編碼RNA (lincRNA)基因是lncRNA中的一類，具有重組染色質和基因組結構、RNA穩定和轉錄調控(包括增強相關活性)等功能。有一些lincRNAs含有小開放閱讀框(smORFs)和編碼功能肽，也被歸類為編碼RNA。lincRNAs可以廣泛地通過多種基因來微調鄰近基因的表達，具有顯著的組織特異性。

第二個問題：如何獲得GTF文件以及如何得到lncRNA對應的GENEID?

解決問題：如下，需要一定的linux基礎知識

2.從GTF文件注釋

想要獲得GTF文件，又需要補充和回顧下下面的知識點

超精華生信ID總結，想踏入生信大門的你-值得擁有

https://github.com/jmzeng1314/bioinformatics123/blob/master/appendix/database.md

https://vip.biotrainee.com/d/165-gtf-gff

基因組注釋文件(GFF,GTF)下載的四種方法

  • [NCBI]
  • [Ensembl]
  • [UCSC]
  • [GeneCode]
（1）不同注釋版本之間的對應關係

hg19，GRCH37和ensembl75是三種國際生物資訊學資料庫資源收集存儲單位，即NCBI，UCSC和ENSEMBL各自發布的人類基因組資訊。 hg系列，hg18/19/38來自UCSC也是目前使用頻率最高的基因組。hg38，是目前的最新版本。 基因組各種版本對應關係綜合來看如下所示 ： NCBI (UCSC) : ENSEMBL

  • GRCh36 (hg18): ENSEMBL release_52.
  • GRCh37 (hg19): ENSEMBL release_59/61/64/68/69/75.
  • GRCh38 (hg38): ENSEMBL release_76/77/78/80/81/82.

ENSEMBL的版本特別複雜也很容易搞混，UCSC的版本就簡單很多，常用的是hg19，最新版本為hg38。

==Ensembl，還有2個大名鼎鼎的計劃就是ENCODE和GENCODE==

Ensembl與ENCODE以及GENCODE計劃之間的關係 Ensembl是ENCODE計劃的子項目。而GENCODE計劃（由Sanger研究所維護）則是ENCODE項目的衍生品，它的目標是為ENCODE項目提供可用的人類基因組和小鼠基因組注釋。Ensembl在ENCODE計劃中的作用是，為人類基因組的組裝提供電腦的自動注釋資訊，並且把這些自動注釋的資訊和來自HAVANA的人工注釋資訊進行合併。GENCODE中的人類和小鼠的基因組注釋和Ensembl資料庫是同步發行的。

GTF和GFF格式是Sanger研究所定義，是一種簡單的、方便的對於DNA、RNA以及蛋白質序列的特徵進行描述的一種數據格式。

==GFF==全稱為general feature format，這種格式主要是用來注釋基因組。

==GTF==全稱為gene transfer format，主要是用來對基因進行注釋。

比如序列的哪裡到哪裡是基因，是轉錄本，是外顯子，內含子或者CDS等等，已經成為序列注釋的通用格式，許多軟體都支援輸入或者輸出gff格式。gff是由tab鍵隔開的9列組成，以下是各列的說明：

  1. seq_id：序列的編號，一般為chr或者contig ID；
  2. source: 注釋的來源，一般為資料庫或者注釋的機構，如果未知，則用點「.」代替
  3. type: 注釋資訊的類型，比如Gene、cDNA、mRNA、CDS等;
  4. start: 該基因或轉錄本在參考序列上的起始位置；
  5. end: 該基因或轉錄本在參考序列上的終止位置；
  6. score: 得分，數字，是注釋資訊可能性的說明，可以是序列相似性比對時的E-values值或者基因預測是的P-values值，.表示為空;
  7. strand: 該基因或轉錄本位於參考序列的正鏈(+)或負鏈(-)上;
  8. phase: 僅對注釋類型為「CDS」有效，表示起始編碼的位置，有效值為0、12. (對於編碼蛋白質的CDS來說，本列指定下一個密碼子開始的位置。每3個核苷酸翻譯一個氨基酸，從0開始，CDS的起始位置，除以3，餘數就是這個值，，表示到達下一個密碼子需要跳過的鹼基個數。該編碼區第一個密碼子的位置，取值0,1,2。0表示該編碼框的第一個密碼子第一個鹼基位於其5』末端；1表示該編碼框的第一個密碼子的第一個鹼基位於該編碼區外；2表示該編碼框的第一個密碼子的第一、二個鹼基位於該編碼區外；如果Feature為CDS時，必須指明具體值；
  9. attributes: 一個包含眾多屬性的列表，格式為「標籤＝值」(tag=value)，以多個鍵值對組成的注釋資訊描述，鍵與值之間用「=」，不同的鍵值用「；」隔開，一個鍵可以有多個值，不同值用「,」分割。如果描述中包括tab鍵以及「，= ；」，要用URL轉義規則進行轉義。
（2）GFF示例==
gff-version 3  description: evidence-based annotation of the human genome (GRCh38), version 31 (Ensembl 97) - long non-coding RNAs  provider: GENCODE  contact: [email protected]  format: gff3  date: 2019-06-27  sequence-region chr1 1 248956422  chr1    HAVANA  gene    29554   31109   .       +       .       ID=ENSG00000243485.5;gene_id=ENSG00000243485.5;gene_type=lncRNA;gene_name=MIR1302-2HG;level=2;hgnc_id=HGNC:52482;tag=ncRNA_host;havana_gene=OTTHUMG0000000095  chr1    HAVANA  transcript      29554   31097   .       +       .       ID=ENST00000473358.1;Parent=ENSG00000243485.5;gene_id=ENSG00000243485.5;transcript_id=ENST00000473358.1;gene_type=lncRNA;gene_name=MIR1302-2HG;transc  chr1    HAVANA  exon    29554   30039   .       +       .       ID=exon:ENST00000473358.1:1;Parent=ENST00000473358.1;gene_id=ENSG00000243485.5;transcript_id=ENST00000473358.1;gene_type=lncRNA;gene_name=MIR1302-2HG;transcr

gencode.v31.long_noncoding_RNAs.gff3.gz

可以看到第9列其實是可以無限擴展的，==只需要以；進行分割即可==。

（3）GTF 與GFF的差異==

GTF文件以及GFF文件都由9列數據組成，這兩種文件的前8列都是相同的（一些小的差別），它們的差別重點在第9列。 GTF文件的第9列同GFF文件不同，雖然同樣是標籤與值配對的情況，但標籤與值之間以空格分開,而不是GFF裡面的=號。相同的是每個特徵之後都要有分號;分割。

（4）GTF示例==
##description: evidence-based annotation of the human genome (GRCh38), version 31 (Ensembl 97) - long non-coding RNAs  ##provider: GENCODE  ##contact: [email protected]  ##format: gtf  ##date: 2019-06-27  chr1    HAVANA  gene    29554   31109   .       +       .       gene_id "ENSG00000243485.5"; gene_type "lncRNA"; gene_name "MIR1302-2HG"; level 2; hgnc_id "HGNC:52482"; tag "ncRNA_host"; havana_gene "OTTHUMG00000000959.2"  chr1    HAVANA  transcript      29554   31097   .       +       .       gene_id "ENSG00000243485.5"; transcript_id "ENST00000473358.1"; gene_type "lncRNA"; gene_name "MIR1302-2HG"; transcript_type "lncRNA"; transcript_nam  chr1    HAVANA  exon    29554   30039   .       +       .       gene_id "ENSG00000243485.5"; transcript_id "ENST00000473358.1"; gene_type "lncRNA"; gene_name "MIR1302-2HG"; transcript_type "lncRNA"; transcript_name "MIR13

gencode.v31.long_noncoding_RNAs.gtf.gz

（5）基本概念的理解

參考一些網路教程，當然推薦大家看書咯。

基於轉錄本上的概念

==外顯子和內含子==： 基因 DNA 分為編碼區和非編碼區，編碼區包含外顯子和內含子，一般非編碼區具有基因表達的調控功能，如啟動子在非編碼區。編碼區則轉錄為 mRNA 並最終翻譯成蛋白質。 外顯子和內含子都被轉錄到 mRNA 前體 hnRNA 中，當 hnRNA 進行剪接變為成熟的 mRNA 時，內含子被切除，而外顯子保留。實際上真正編碼蛋白質的是外顯子，而內含子則無編碼功能。 內含子存在於DNA 中，在轉錄的過程中，DNA 上的內含子也會被轉錄到前體 RNA 中，但前體 RNA 上的內含子會在 RNA 離開細胞核進行翻譯前被切除。 promoter不屬於intron和exon的任何一個，屬於noncoding sequence. ==開放讀碼框====ORF== 開放讀碼框是從一個起始密碼子開始到一個終止密碼子結束的一段序列；不是所有讀碼框都能被表達出蛋白產物，或者能表達出佔有優勢或者能產生生物學功能的蛋白。內含子和外顯子指的就是一個開放閱讀框（ORF）內編碼的部分和不編碼的部分。

基於翻譯上的概念

mRNA 包括 UTR 和 CDS ==UTR==：UTR,Untranslated Regions一班指的是一個轉錄本（transcript）3'和5'不參與編碼的區域即非翻譯區，是信使RNA（mRNA）分子兩端的非編碼片段。UTR區不參與編碼，但是不是說他們沒有功能，只是不被翻譯成具有功能的蛋白質。多數基因都有UTR,它們也是外顯子拼接的產物。UTR在DNA序列中是算外顯子的區域。 ==CDS==：CDS,Sequencecodingfor aminoacids in protein 蛋白質編碼區 ，CDS 是 Coding sequence的縮寫，是編碼一段蛋白產物的序列，CDS 必定是一個 ORF 。但也可能包括很多 ORF 。反之，每個 ORF 不一定都是 CDS 。 外顯子與 CDS 區不是完全一致的，CDS 區一定屬於外顯子，但是外顯子不一定是 CDS 區，也就是說外顯子不一定都能翻譯成蛋白。

三張圖片幫助理解

外顯子和內含子是轉錄本上面的概念，是基於轉錄這個行為的定義。 而5端UTR區域和CDS區域，還有3端UTR區域,是基於翻譯這種行為的定義！

（6）如何獲得GTF、GFF文件

上面是關於基本概念的理解，那麼從哪裡下載GFF文件和GTF文件？

這邊是從gencode官網，步驟如下：

登錄gencode官網：https://www.gencodegenes.org/

1.點擊How to access data

2.點擊human

3.選擇最新的版本version 31

4.下載GTF或GFF文件

接下來就可以解決上面的第二個問題，即如何獲得GTF文件以及如何得到lncRNA對應的GENEID?

（7）獲得lncRNA對應的GENEID

下載的gencode.v31.long_noncoding_RNAs.gtf.gz文件用FileZilla上傳到伺服器，輸入如下

zless -S gencode.v31.long_noncoding_RNAs.gtf.gz|grep -w 'gene'|cut -f 9|cut -d ';' -f 1-3|tr ';' ' '|awk '{print $2 $4 $6}'|less -S > anno.txt

再將anno.txt用FileZilla傳到本地。

（8）在R里讀取id注釋文件
######step3-id_convertion  options(stringsAsFactors = F)  lncRNA_genecode<-read.table('anno.txt')  tmp<-sapply(strsplit(lncRNA_genecode$V1,'.',fixed = T),function(x)x[1])  lncRNA_genecode$GENEID<-tmp  colnames(lncRNA_genecode)[2]<-'GENEBIOTYPE'  colnames(lncRNA_genecode)[3]<-'gene_name'  table(lncRNA_genecode$GENEID %in% gene2sym$GENEID)  tmp1<-lncRNA_genecode[lncRNA_genecode$GENEID %in% gene2sym$GENEID,]  lncRNA<-gene2sym[match(tmp1$GENEID,gene2sym$GENEID),]    tmp<-sapply(strsplit(rownames(dat),'.',fixed = T),function(x)x[1])  rownames(dat)<-tmp  mRNA<-gene2sym[which(gene2sym$GENEBIOTYPE=="protein_coding"),]  dim(mRNA)

最終得到的mRNA和lncRNA的GENEID(也就是ENSG開頭的)和SYMBOL的對應關係矩陣就得到了，如下圖

最後，得到上面的mRNA和lncRNA的GENEID和symbol的對應關係，在通過GENEID對應到表達矩陣——step1得到的dat，分別走標準差異分析流程，就可以分別獲得mRNA和lncRNA中的差異基因啦！這樣就可以很容易得到文章中的第一張圖片啦

後記

可以看到，僅僅是復現文章的第一章圖表就已經非常的不容易了，說起來僅僅是把基因分類，實際上要理解背後的生物學知識，資料庫來源知識，需要linux和R能力，以前我覺得學習很簡單，但是看到了太多的零基礎和負基礎，負起點的朋友的掙扎才發現完全不是那麼一回事。

當然了，共資料庫挖掘需要的基礎linux和r技巧好好掌握。

Tags:

VirMach 便宜 VPS