GTF文件有什么用啊？别的不谈，最起码能提lncRNA

[ixxmu]

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

[ixxmu]

GTF文件有什么用啊？别的不谈，最起码能提lncRNA by 果子学生信

我们直接去下载他， http://asia.ensembl.org/index.html

解压缩后，使用R语言可以读取，有多中方案，最终我选择的是rtracklayer::import，而且我认为这是最好的方法。

首先我们来安装rtracklayer这个包

if(length(getOption("CRAN"))==0) options(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")
if(!require("rtracklayer")) BiocManager::install("rtracklayer")

直接读取，大概要2分钟，然后转换成data.frame格式

gtf1 <- rtracklayer::import('Homo_sapiens.GRCh38.94.chr.gtf')
gtf_df <- as.data.frame(gtf1)

这是一行270万行，27列的文件,可能是普通人接触到的最大数据，里面记录了基因的构成信息，每个基因有在基因组上的位置，有多少转录本，是否编码。

dim(gtf_df)
[1] 2736845      27

我们把这个数据保存成Rdata格式，这就是课堂上抽取lncRNA，mRNA以及名称转换需要的gtf文件，使用的时候load即可

save(gtf_df,file = "gtf_df.Rda")

我们通过行名来了解一下他,有27个。

colnames <- data.frame(colnames=colnames(gtf_df))
colnames
1 seqnames
2 start
3 end
4 width
5 strand
6 source
7 type
8 score
9 phase
10 gene_id
11 gene_version
12 gene_name
13 gene_source
14 gene_biotype
15 transcript_id
16 transcript_version
17 transcript_name
18 transcript_source
19 transcript_biotype
20 tag
21 transcript_support_level
22 exon_number
23 exon_id
24 exon_version
25 protein_id
26 protein_version
27 ccds_id

然后我们就可以利用这个文件结合读取好的TCGA counts文件，抽离出mRNA和lncRNA 其中exprdfnopoint就是去掉ensemble id后面小数点的表达文件，ensemble id是第一列，每个样本是一列。大概是这样的：

获取mRNA

library(dplyr)
library(tidyr)
## mRNA
mRNA_exprSet <- gtf_df %>% 
  dplyr::filter(type=="gene",gene_biotype=="protein_coding") %>% #筛选gene,和编码指标
  dplyr::select(c(gene_name,gene_id,gene_biotype)) %>% 
  dplyr::inner_join(expr_df_nopoint,by ="gene_id") %>% 
  tidyr::unite(gene_id,gene_name,gene_id,gene_biotype,sep = " | ")

获取lncRNA

library(dplyr)
library(tidyr)
# lncRNA
ncRNA <- c("sense_overlapping","lincRNA","3prime_overlapping_ncRNA",
           "processed_transcript","sense_intronic",
           "bidirectional_promoter_lncRNA","non_coding",
           "antisense_RNA")
LncRNA_exprSet <- gtf_df %>% 
  dplyr::filter(type=="transcript",transcript_biotype %in% ncRNA) %>% #注意这里是transcript_biotype
  dplyr::select(c(gene_name,gene_id,transcript_biotype)) %>% 
  dplyr::distinct() %>% #删除多余行
  dplyr::inner_join(expr_df_nopoint,by ="gene_id") %>% 
  tidyr::unite(gene_id,gene_name,gene_id,transcript_biotype,sep = " | ")

除此以外，还能干什么呢？ 比如我想知道21号染色体上有多少个基因？可以的

## 21号染色体上的编码基因
gene21 <- gtf_df %>% 
  dplyr::filter(seqnames=="21",type=="gene",gene_biotype=="protein_coding") %>%
  dplyr::select(c(gene_name,gene_id)) 

我们发现有231个编码基因

人类基因组上每个染色体上有多少个基因，知道么？可以的。

## 每个染色体上有多少个基因
gene_chr <- gtf_df %>% 
  dplyr::filter(type=="gene")
data <- table(gene_chr$seqnames)
barplot(data,col = data)

好了，明天见。