转录组测序后的表达量矩阵是否需要清理非编码基因呢

生信技能树

发布于 2024-11-21 01:54:15

13100

代码可运行

运行总次数：0

代码可运行

在前面的笔记里面：有一些错误在图片上面显示不出来，我们提到了一个数据挖掘文章把tcga数据库的肝癌转录组测序数据集里面的差异分析弄反了，方法学描述是：using the package DEG- seq2, Adj. p value < 0.05 and |logFC| > 2 were regarded as the cut-off criteria. This identified 2162 genes met the standards

实际上我们很容易去复现，但是很难拿到同样的差异情况。

首先看看表达量矩阵

# 魔幻操作，一键清空
rm(list = ls()) 
options(stringsAsFactors = F)
library(data.table)
a1=fread( 'input/TCGA-LIHC.htseq_counts.tsv.gz' , 
          data.table = F)
dim(a1)
a1[1:4,1:4]
a1[(nrow(a1)-5):nrow(a1),1:4]
dim(a1)
# all data is then log2(x+1) transformed. 
#length(unique(a1$AccID))
#length(unique(a1$GeneName))
mat= a1[,2:ncol(a1)] 
mat[1:4,1:4] 
mat=mat[1:(nrow(a1)-4),]
mat=ceiling(2^(mat)-1) #log2(x+1) transformed.
mat[1:4,1:4] 
rownames(mat) = gsub('[.][0-9]+','',a1$Ensembl_ID[1:(nrow(a1)-4)])
keep_feature <- rowSums (mat > 1) > 1
colSums(mat)/1000000
table(keep_feature)
mat <- mat[keep_feature, ]
mat[1:4,1:4] 
mat=mat[,  colSums(mat)/1000000 >10]
dim(mat) 
colnames(mat) 
ensembl_matrix=mat
colnames(ensembl_matrix) 
ensembl_matrix[1:4,1:4]

library(AnnoProbe)
head(rownames(ensembl_matrix))
ids=annoGene(rownames(ensembl_matrix),'ENSEMBL','human')
head(ids)
sort(table(ids$biotypes))

可以看到，表达量矩阵里面的基因分类非常复杂：

> tail(as.data.frame(sort(table(ids$biotypes))))
                     Var1  Freq
31                    TEC   829
32               misc_RNA   978
33 unprocessed_pseudogene  1395
34   processed_pseudogene  7189
35                 lncRNA 12009
36         protein_coding 19008

但是蛋白质编码基因就不到2万个，数量比较多的还有：

**Unprocessed Pseudogene (未处理的伪基因)**：

伪基因（Pseudogene）是基因组中的一种非功能性基因，它们通常来源于功能性基因的复制，但在进化过程中失去了原有的功能。
未处理的伪基因是指那些保留了与原始功能性基因相似的序列，但缺少正常的剪接位点，因此不能产生成熟的mRNA。

**Processed Pseudogene (已处理的伪基因)**：

已处理的伪基因是指那些经历了类似功能性基因的剪接过程，但最终不产生功能性蛋白质的基因。
它们通常通过逆转录事件产生，其中功能性基因的mRNA被逆转录并插入基因组中，但缺乏启动子和调控元件，因此不能被有效转录。

**Long Non-Coding RNA (长链非编码RNA, lncRNA)**：

lncRNA是一类长度超过200个核苷酸的RNA分子，它们不编码蛋白质，但在细胞中扮演多种调控角色。
lncRNA可以参与基因表达的调控，如染色质重塑、转录激活或抑制、RNA剪接和翻译调控等。

**Protein Coding (蛋白质编码)**：

蛋白质编码基因是指那些能够转录并翻译成蛋白质的基因。
这些基因的转录产物（mRNA）经过剪接和翻译，最终形成具有特定功能的蛋白质。

而且很明显的可以看到，两万个蛋白质编码基因其实是占据了每个样品的90%的表达量，其它三四万个基因反而是表达量微乎其微！

  library(AnnoProbe)
  gs=annoGene(rownames(ensembl_matrix),'ENSEMBL','human')
  head(gs)
  tail(sort(table(gs$biotypes)))
  pd_genes=gs[gs$biotypes=='protein_coding',3]
  non_genes=gs[gs$biotypes !='protein_coding',3]
  pd_matrix=ensembl_matrix[rownames(ensembl_matrix) %in% pd_genes,]
  non_matrix=ensembl_matrix[rownames(ensembl_matrix) %in% non_genes,] 
  
  
>   as.numeric(fivenum(colSums(pd_matrix)/1e6))
[1]  19.81015  43.34670  50.60327  57.11173 118.90806
>   as.numeric(fivenum(colSums(non_matrix)/1e6))
[1] 0.3917600 0.8467285 1.0720895 1.3846255 9.5381310

如果是使用全部的基因走后面的转录组测序矩阵差异分析，可以看到上调基因会偏多：

上调基因会偏多

如果是过滤一下，上调基因可以砍掉一大半，而下调基因影响才三分之一：

上调基因可以砍掉一大半，而下调基因影响才三分之一

对比两次差异分析结果

首先呢，几乎是不会影响两万多个蛋白质编码基因的差异分析结果，因为两次差异分析的变化情况是几乎是一模一样：

rm(list = ls())
library(data.table)
load('../deg-all-genes/DEG_deseq2.Rdata')  
DEG_se=DEG_deseq2
load('../deg-for-protein/DEG_deseq2.Rdata') 
DEG_pe=DEG_deseq2
 
ids=intersect(rownames(DEG_se),
              rownames(DEG_pe))
df= data.frame(
  DEG_se = DEG_se[ids,'log2FoldChange'],
  DEG_pe = DEG_pe[ids,'log2FoldChange']
)
library(ggpubr) 
ggscatter(df, x = "DEG_se", y = "DEG_pe",
          color = "black", shape = 21, size = 1, # Points color, shape and size
          add = "reg.line",  # Add regressin line
          add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
          conf.int = TRUE, # Add confidence interval
          cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
          cor.coeff.args = list(method = "pearson",  label.sep = "\n")
)

而且呢，虽然上下调基因数量不一样，但是也很难说有什么冲突：

> DEG_se=modify_deg(DEG_se) 
> table(DEG_se$regulated)

  down normal     up 
  1865  23980   6983 
> DEG_pe=modify_deg(DEG_pe) 
> table(DEG_pe$regulated)

  down normal     up 
  1232  14477   3278 
> ids=intersect(rownames(DEG_se),
+               rownames(DEG_pe))
> df= data.frame(
+   DEG_se = DEG_se[ids,'regulated'],
+   DEG_pe = DEG_pe[ids,'regulated']
+ )
> table(df)
        DEG_pe
DEG_se    down normal    up
  down    1224     28     0
  normal     0  13430    65
  up         0      2  3213
>