文章/答案/技术大牛

发布

社区首页 >专栏 >就想把表达矩阵区分成为蛋白编码基因和非编码有这么难吗？

就想把表达矩阵区分成为蛋白编码基因和非编码有这么难吗？

生信技能树

发布于 2019-09-24 06:01:25

4K00

代码可运行

文章被收录于专栏：生信技能树生信技能树

运行总次数：0

代码可运行

大家好，我是梦梦，以为自己半年前完成了数据挖掘培训学习，结果最近想看看自己是否够格成为优秀学徒，才发现自己与他们差距有多大，所以把自己思考过程分享一下，希望大家也都实践起来，不要自认为懂了就懂了！前面生信技能树公众号优秀学徒发布了七步走纯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列组成，以下是各列的说明：

seq_id：序列的编号，一般为chr或者contig ID；
source: 注释的来源，一般为数据库或者注释的机构，如果未知，则用点“.”代替
type: 注释信息的类型，比如Gene、cDNA、mRNA、CDS等;
start: 该基因或转录本在参考序列上的起始位置；
end: 该基因或转录本在参考序列上的终止位置；
score: 得分，数字，是注释信息可能性的说明，可以是序列相似性比对时的E-values值或者基因预测是的P-values值，.表示为空;
strand: 该基因或转录本位于参考序列的正链(+)或负链(-)上;
phase: 仅对注释类型为“CDS”有效，表示起始编码的位置，有效值为0、12. (对于编码蛋白质的CDS来说，本列指定下一个密码子开始的位置。每3个核苷酸翻译一个氨基酸，从0开始，CDS的起始位置，除以3，余数就是这个值，，表示到达下一个密码子需要跳过的碱基个数。该编码区第一个密码子的位置，取值0,1,2。0表示该编码框的第一个密码子第一个碱基位于其5’末端；1表示该编码框的第一个密码子的第一个碱基位于该编码区外；2表示该编码框的第一个密码子的第一、二个碱基位于该编码区外；如果Feature为CDS时，必须指明具体值；
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: gencode-help@ebi.ac.uk
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: gencode-help@ebi.ac.uk
##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)