2023-01-09 用GenomicFeatures包計算基因長度

本文包括內(nèi)容：

1. 基因長度定義
2. 常用注釋包簡介
3. 基因長度計算:用GenomicFeatures包 計算TxDb.Hsapiens.UCSC.hg38.knownGene中四種不同定義的基因長度。
參考： 
https://blog.csdn.net/Candle_light/article/details/90598308
https://mp.weixin.qq.com/s?__biz=MzAxMDkxODM1Ng==&mid=2247488656&idx=1&sn=c2472cca80a0b0f8b831fbf5c69f2a34&chksm=9b48542bac3fdd3dd99474c42b2af092e462842007e21b28724936a3c9fc9f9c14625d355fd1&scene=21#wechat_redirect
https://mp.weixin.qq.com/s?__biz=MzU4NjU4ODQ2MQ==&mid=2247490627&idx=1&sn=f23242af5baa6cd6c07ff3558d65c97b&chksm=fdf85401ca8fdd17017107413706c4c259c91385ddc17e0778539d327f2f1d49f0f5072ea688&scene=21#wechat_redirect

1. 基因長度定義

非冗余外顯子(EXON)長度之和;
挑選基因的最長轉(zhuǎn)錄本 ;
選取多個轉(zhuǎn)錄本長度的平均值;
非冗余 CDS 長度之和

2. 常用注釋包簡介

2.1 TxDb包

包含位置信息。包里的注釋內(nèi)容，可以通過包名來判斷。包名的格式是：TxDb.Species.Source.Build.Table
例如，TxDb.Hsapiens.UCSC.hg19.knownGene，代表：
物種是 Homo sapiens
來源于 UCSC genome browser
基因組版本是 hg19 (their version of GRCh37)
包含的注釋內(nèi)容是 knownGene table

1.2 EnsDb包

EnsDb和TxDb類似，只不過它是基于Ensembl數(shù)據(jù)庫的
例如：EnsDb.Hsapiens.v79

3. GRanges，Ranges

Ranges 對象功能十分強大，可以讓我們基于基因組位置信息輕松選取出相關(guān)數(shù)據(jù)。GRanges和GRangesLists對象，就類似于data.frame和List數(shù)據(jù)結(jié)構(gòu)，我們可以使用[符號去選取子集。

3.1 GRanges

注釋包比較常用的一個功能，是將注釋包中的位置信息提取出來，添加到GRanges或者GRangesList對象中。

  • 將GRanges對象寫成txt
exon_txdb=exons(txdb)
tmp=as.data.frame(exon_txdb)
write.table(tmp,"exon.pos",row.names=F)

genes_txdb=genes(txdb)
tmp=as.data.frame(genes_txdb)
write.table(tmp,"gene.pos",row.names=F)
  • 其他操作
seqnames(exon_txdb) 
ranges(exon_txdb) #返回外顯子的起始終止位點，長度，以及其它信息
strand(exon_txdb)#返回外顯子的正負鏈信息，要么在正鏈要么在負鏈
mcols(exon_txdb)#返回exon的id編號，1到724366個
seqlengths(exon_txdb)#返回每條染色體的長度信息 names,length
GRanges對象還有很多其它類型的操作，非常好玩的，split,shift,resize,flank,reduce,gaps,disjoin,coverage
其它求交集并集和都可以用,union,intersect,setdiff,pintersect,psetdiff

3.2 GRangesLists

以List形式呈現(xiàn)，每個基因所對應(yīng)的轉(zhuǎn)錄本，在基因組上的坐標(biāo)

3.3 其他

transcripts(), genes()函數(shù)選取位置信息， 編碼區(qū)可以使用cds(),promoters(),exons()來選取

library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene 
genes(txdb)#29794個基因
exons(txdb)#724366 個exon
transcripts(txdb)  #272352個轉(zhuǎn)錄本
cds(txdb)#306522個編碼區(qū)
transcriptsBy(txdb,by="gene")#對象按照gene來對轉(zhuǎn)錄本分組
exonsBy,cdsBy來對它進行處理，都會返回Granges對象

4. 根據(jù)四種定義計算基因長度

目前了解到兩種方法:1) 使用R包GenomicFeatures 2)使用R包 rtracklayer；

首先，載入txdb文件

##創(chuàng)建txdb文件
if(!require(GenomicFeatures))BiocManager::install("GenomicFeatures")
  library(GenomicFeatures)
  txdb <- makeTxDbFromGFF("gencode.v36.annotation.gtf.gz",format="gtf")
##或者下載
BiocManager::install("TxDb.Hsapiens.UCSC.hg38.knownGene")
library(TxDb.Hsapiens.UCSC.hg38.knownGene) 
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene 

#### R包GenomicFeatures####
library(GenomicFeatures)
library(TxDb.Hsapiens.UCSC.hg38.knownGene) 
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene 

##① 非冗余外顯子(EXON)長度之和 
exons.list.per.gene <- exonsBy(txdb,by="gene")
head(as.data.frame(exons.list.per.gene))
  group group_name seqnames    start      end width strand exon_id
1     1          1    chr19 58345178 58347029  1852      -  603914
2     1          1    chr19 58345183 58347029  1847      -  603915
3     1          1    chr19 58346854 58347029   176      -  603916
4     1          1    chr19 58346858 58347029   172      -  603917
5     1          1    chr19 58346860 58347029   170      -  603918
6     1          1    chr19 58347353 58347634   282      -  603919
  exon_name
1      <NA>
2      <NA>
3      <NA>
4      <NA>
5      <NA>
6      <NA>
# 然后計算每個基因的非冗余外顯子長度(widths)并加和。
exonic.gene.sizes <- sum(width(GenomicRanges::reduce(exons.list.per.gene)))
          exonic.gene.sizes
head(as.data.frame(exonic.gene.sizes))
1                      6528
10                     1834
100                   12177
1000                   8583
10000                 22459
100008586               693
gle <- data.frame(gene_id=names(exonic.gene.sizes),
                    length=exonic.gene.sizes) ##31456個基因
            gene_id length
1                 1   6528
10               10   1834
100             100  12177
1000           1000   8583
10000         10000  22459
100008586 100008586    693
##②非冗余cds長度之和
cds.list.per.gene <- cdsBy(txdb,by="gene")
head(as.data.frame(cds.list.per.gene))
  group group_name seqnames    start      end width strand cds_id
1     1          1    chr19 58347022 58347029     8      - 252822
2     1          1    chr19 58347353 58347640   288      - 252823
3     1          1    chr19 58350370 58350651   282      - 252824
4     1          1    chr19 58350594 58350651    58      - 252825
5     1          1    chr19 58351391 58351687   297      - 252826
6     1          1    chr19 58352098 58352184    87      - 252827
  cds_name
1     <NA>
2     <NA>
3     <NA>
4     <NA>
5     <NA>
6     <NA>

# 然后計算每個基因的非冗余外顯子長度(widths)并加和。
cds.gene.sizes <- sum(width(GenomicRanges::reduce(cds.list.per.gene)))
head(as.data.frame(cds.gene.sizes))
          cds.gene.sizes
1                   1575
10                   873
100                 1979
1000                2978
10000               3296
100008586            354
glc <- data.frame(gene_id=names(cds.gene.sizes),
                  length=cds.gene.sizes)#19626個基因
            gene_id length
1                 1   1575
10               10    873
100             100   1979
1000           1000   2978
10000         10000   3296
100008586 100008586    354
##③挑選基因的最長轉(zhuǎn)錄本
t_l=transcriptLengths(txdb) #272352
t_l=na.omit(t_l)#229531
t_l=t_l[order(t_l$tx_len,decreasing = T),]# 這里把同樣的基因的多個轉(zhuǎn)錄本長度排序了
head(t_l)
        tx_id            tx_name gene_id nexon tx_len
38851   38851  ENST00000589042.5    7273   363 109224
38852   38852  ENST00000591111.5    7273   313 104301
38849   38849 ENST00000342992.11    7273   312 101520
138913 138913  ENST00000597346.1   10984     1  91667
38850   38850  ENST00000460472.6    7273   191  82029
38853   38853 ENST00000342175.11    7273   191  81357
t_l=t_l[!duplicated(t_l$gene_id),]#去重，直接保留轉(zhuǎn)錄本最長的
t_l=t_l[order(t_l$gene_id,decreasing = F),]
glm=data.frame(gene_id=t_l$gene_id,length=t_l$tx_len)
head(glm)
    gene_id length
1         1   3382
2        10   1285
3       100   5965
4      1000   4016
5     10000   7950
6 100008586    585
##④ 選取多個轉(zhuǎn)錄本長度的平均值 
t_l=transcriptLengths(txdb) #272352
t_l=na.omit(t_l)#229531
t_l=t_l[order(t_l$gene_id,decreasing = T),]
gla <- t_l %>% group_by(gene_id) %>% 
summarise(tx_ave_len = round(mean(tx_len),0))#31456個基因
gla <-as.data.frame(gla)
head(gla)
    gene_id tx_ave_len
1         1       1842
2        10        850
3       100       1708
4      1000       2744
5     10000       2436
6 100008586        573

比較不同方法得到的基因長度:

a = gle %>% inner_join(glc,"gene_id") %>% 
  inner_join(glm,"gene_id") %>% 
  inner_join(gla,"gene_id")
colnames(a) = c("gene_id","exon","CDS","max","average") # 19626個基因
head(a)
    gene_id  exon  CDS  max average
1         1  6528 1575 3382    1842
2        10  1834  873 1285     850
3       100 12177 1979 5965    1708
4      1000  8583 2978 4016    2744
5     10000 22459 3296 7950    2436
6 100008586   693  354  585     573
boxplot(log2(a[,-1]),outline = F)
不同計算方法基因長度

gene_id 轉(zhuǎn)換到symbol看下

library(org.Hs.eg.db)
s2g=toTable(org.Hs.egSYMBOL)
head(s2g)
aa =merge(s2g,a,by='gene_id') #19601個基因，比id轉(zhuǎn)換前少了25個基因
head(aa)
    gene_id  symbol  exon  CDS  max average
1         1    A1BG  6528 1575 3382    1842
2        10    NAT2  1834  873 1285     850
3       100     ADA 12177 1979 5965    1708
4      1000    CDH2  8583 2978 4016    2744
5     10000    AKT3 22459 3296 7950    2436
6 100008586 GAGE12F   693  354  585     573

1.不同版本的注釋文件，得到的長度不太一樣，我的結(jié)果是來自TxDb.Hsapiens.UCSC.hg38.knownGene，小潔的結(jié)果是Homo_sapiens.GRCh38.103.chr.gtf.gz。 可以看到exon 和cds 的計算方式差別更小，最長轉(zhuǎn)錄本和平均轉(zhuǎn)錄本的計算方式差的比較多。


我的結(jié)果

小潔的結(jié)果

2. boxplot的結(jié)果和小潔的也不太一樣？ 我的看起來exon計算結(jié)果基因長度最長，小潔的是按最長轉(zhuǎn)錄本的方式，基因長度最長.

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容