RNA-seq入门实战（三）：在R里面整理表达量counts矩阵

连续两次求贤令：曾经我给你带来了十万用户，但现在祝你倒闭，以及 生信技能树知识整理实习生招募，让我走大运结识了几位优秀小伙伴！大家开始根据我的ngs组学视频进行一系列公共数据集分析实战，其中几个小伙伴让我非常惊喜，不需要怎么沟通和指导，就默默的完成了一个实战！

他前面的分享是：

下面是他对我们b站转录组视频课程的详细笔记

本节概览：

1. 从featureCounts输出文件中获取counts与TPM矩阵:
   读取counts.txt构建counts矩阵；样品的重命名和分组；counts与TPM转换；基因ID转换；初步过滤低表达基因与保存counts数据

2. 从salmon输出文件中获取counts与TPM矩阵:
   用tximport包读取quant.sf构建counts与TPM矩阵；样品的重命名和分组；初步过滤低表达基因与保存counts数据

承接上节RNA-seq入门实战（二）：上游数据的比对计数——Hisat2与Salmon之前已经得到了featureCounts与Salmon输出文件（counts、salmon）和基因ID转化文件（g2s_vm25_gencode.txt、t2s_vm25_gencode.txt)。

一般为了对样品进行分组注释我们还需要在GEO网站下载样品Metadata信息表SraRunTable.txt，接下来就需要在R中对输出结果进行操作，转化为我们想要的基因表达counts矩阵。

image.png

一、从featureCounts输出文件中获取counts矩阵

1. 读取counts.txt构建counts矩阵，进行样品的重命名和分组

###环境设置
rm(list=ls())
options(stringsAsFactors = F) 
library(tidyverse) # ggplot2 stringer dplyr tidyr readr purrr  tibble forcats
library(data.table) #多核读取文件
setwd("C:/Users/Lenovo/Desktop/test")

#### 对counts进行处理筛选得到表达矩阵 ####
a1 <- fread('./counts/counts.txt',
              header = T,data.table = F)#载入counts，第一列设置为列名
colnames(a1)
counts <- a1[,7:ncol(a1)] #截取样本基因表达量的counts部分作为counts 
rownames(counts) <- a1$Geneid #将基因名作为行名
#更改样品名
colnames(counts)
colnames(counts) <- gsub('/home/test/align/bam/','', #删除样品名前缀
                   gsub('_sorted.bam','',  colnames(counts))) #删除样品名后缀


#### 导入或构建样本信息,  进行列样品名的重命名和分组####
b <- read.csv('./SraRunTable.txt')
b
name_list <- b$source_name[match(colnames(counts),b$Run)]; name_list  #选择所需要的样品信息列
nlgl <- data.frame(row.names=colnames(counts),
                   name_list=name_list,
                   group_list=name_list)
fix(nlgl)  #手动编辑构建样品名和分组信息
name_list <- nlgl$name_list
colnames(counts) <- name_list #更改样品名
group_list <- nlgl$group_list
gl <- data.frame(row.names=colnames(counts), #构建样品名与分组对应的数据框
                 group_list=group_list)

这里给样品名加上_1、_2表示重复样品，根据这两类细胞的多能性状态将其分组为naive和primed

fix(nlgl)编辑构建样品名和分组信息

2. counts与TPM转换

基因表达量一般以TPM或FPKM为单位来展示，所以还需要进行，若还想转化为FPKM或CPM可参见Counts FPKM RPKM TPM 的转化与 获取基因有效长度的N种方法

#### counts，TPM转化 ####
# 注意需要转化的是未经筛选的counts原始矩阵
### 从featurecounts 原始输出文件counts.txt中提取Geneid、Length(转录本长度)，计算tpm
geneid_efflen <- subset(a1,select = c("Geneid","Length"))
colnames(geneid_efflen) <- c("geneid","efflen")  

### 取出counts中geneid对应的efflen
efflen <- geneid_efflen[match(rownames(counts),
                              geneid_efflen$geneid),
                        "efflen"]

### 计算 TPM 公式
 #TPM (Transcripts Per Kilobase Million)  每千个碱基的转录每百万映射读取的Transcripts
  counts2TPM <- function(count=count, efflength=efflen){
    RPK <- count/(efflength/1000)   #每千碱基reads (Reads Per Kilobase) 长度标准化
    PMSC_rpk <- sum(RPK)/1e6        #RPK的每百万缩放因子 (“per million” scaling factor ) 深度标准化
    RPK/PMSC_rpk              
  }  

tpm <- as.data.frame(apply(counts,2,counts2TPM))
colSums(tpm)

3. 基因ID转换

若上游中采用的是UCSC的基因组和gtf注释文件，则表达矩阵行名就是我们常见的gene symbol基因名；若上游采用的是gencode或ensembl基因组和gtf注释文件，那么我们就需要将基因表达矩阵行名的Ensembl_id(gene_id)转换为gene symbol (gene_name)了。
在转换时经常会出现多个Ensembl_id对应一个gene symbol的情形，此时就出现了重复的gene symbol。此时就需要我们在进行基因ID转换前去除重复的gene symbol。
下面展示转化ID并合并所有重复symbol的方法，其他基因名去重复方法参见Ensembl_id转换与gene symbol基因名去重复的两种方法 - 简书 (jianshu.com)

#合并所有重复symbol
g2s <- fread('g2s_vm25_gencode.txt',header = F,data.table = F) #载入从gencode的gtf文件中提取的信息文件
colnames(g2s) <- c("geneid","symbol")
  
symbol <- g2s[match(rownames(counts),g2s$geneid),"symbol"] #匹配counts行名对应的symbol
table(duplicated(symbol))  #统计重复基因名
  
###使用aggregate根据symbol列中的相同基因进行合并 
counts <- aggregate(counts, by=list(symbol), FUN=sum)
counts <- column_to_rownames(counts,'Group.1')
  
tpm <- aggregate(tpm, by=list(symbol), FUN=sum) ###使用aggregat 将symbol列中的相同基因进行合并 
tpm <- column_to_rownames(tpm,'Group.1')

id转换前

id转换后

4. 初步过滤低表达基因与保存counts数据

我们的数据中会有很多低表达甚至不表达的基因，在后续分析中可能会影响数据的分析判断，因此需要对低表达的基因进行筛除处理。筛选标准不唯一，依自己数据情况而定。在这里展示筛选出至少在重复样本数量内的表达量counts大于1的行（基因），可以看到超过一半以上的基因都被筛掉了。（这个是正常现象，因为我们的gtf文件里面的基因数量太多了，都是五六万个，而正常情况下我们的样品里面就两万多个基因是有表达量的）

#### 初步过滤低表达基因 ####（筛选标准不唯一、依情况而定）
#筛选出至少在重复样本数量内的表达量counts大于1的行（基因）
keep_feature <- rowSums(counts>1) >= 2
table(keep_feature)  #查看筛选情况，FALSE为低表达基因数（行数），TURE为要保留基因数
#FALSE  TRUE 
#35153 20339 

counts_filt <- counts[keep_feature, ] #替换counts为筛选后的基因矩阵（保留较高表达量的基因）
tpm_filt <- tpm[keep_feature, ]

#### 保存数据 ####
counts_raw=counts #这里重新命名方便后续分析调用
counts=counts_filt
tpm=tpm_filt

save(counts_raw, counts, tpm, 
     group_list, gl,
     file='./1.counts.Rdata')

二、从salmon输出文件中获取counts矩阵

需要用到tximport包从salmon输出文件中获取counts矩阵，在tximport函数中输入quant.sf文件路径、转换类型type = "salmon"、以及转录本与基因名gene symbol对应关系文件（即我们之前得到的t2s_vm25_gencode.txt）就可以转换得到各基因的定量关系了。其他步骤与操作featureCounts输出文件类似。

这里只展示了获取基因表达的TPM值，如果还想了解如何获得FPKM值请参考文章：获取基因有效长度的N种方法中第二部分内容以及Counts FPKM RPKM TPM 的转化。

rm(list=ls())
options(stringsAsFactors = F)
library(tximport) #Import transcript-level abundances and counts for transcript- and gene-level analysis packages
library(tidyverse) # ggplot2 stringer dplyr tidyr readr purrr  tibble forcats
library(data.table) #多核读取文件
setwd("C:/Users/Lenovo/Desktop/test/")


####  salmon原始文件处理  ####
##载入transcript_id和symbol的对应关系文件
t2s <- fread("t2s_vm25_gencode.txt", data.table = F, header = F); head(t2s)

##找到所有quant.sf文件所在路径  导入salmon文件处理汇总
files <- list.files(pattern="*quant.sf",recursive=T, full.names = T); files  #显示目录下所有符合要求的文件
txi <- tximport(files, type = "salmon", tx2gene = t2s)

##提取文件夹中的样品名作为counts行名
cn <- sapply(strsplit(files,'\\/'), function(x) x[length(x)-1]); cn
colnames(txi$counts) <- gsub('_quant','',cn); colnames(txi$counts)

##提取出counts/tpm表达矩阵
counts <- as.data.frame(apply(txi$counts,2,as.integer)) #将counts数取整
rownames(counts) <- rownames(txi$counts) 
tpm <- as.data.frame(txi$abundance)  ###abundance为基因的Tpm值
colnames(tpm) <- colnames(txi$counts)
 

#### 导入或构建样本信息,  进行列重命名和分组 ####
b <- read.csv('./SraRunTable.txt')
b
name_list <- b$source_name[match(colnames(counts),b$Run)]; name_list
nlgl <- data.frame(row.names=colnames(counts),
                 name_list=name_list,
                 group_list=name_list)
fix(nlgl)
name_list <- nlgl$name_list
colnames(counts) <- name_list
colnames(tpm) <- name_list

group_list <- nlgl$group_list
gl <- data.frame(row.names=colnames(counts),
                 group_list=group_list)


#### 初步过滤低表达基因 ####
#筛选出至少在重复样本数量内的表达量counts大于1的行（基因）
keep_feature <- rowSums(counts > 1) >= 2               #ncol(counts)/length(table(group_list)) 
table(keep_feature)  #查看筛选情况
counts_filt <- counts[keep_feature, ] #替换counts为筛选后的基因矩阵（保留较高表达量的基因）
tpm_filt <- tpm[keep_feature, ]


#### 保存数据 ####
counts_raw=counts  
counts=counts_filt
tpm=tpm_filt

save(counts_raw, counts, tpm,
     group_list, gl, txi, #注意保存txi文件用于DESeq2分析
     file='salmon/1.counts.Rdata')

---

通过以上步骤，成功从featureCounts或Salmon输出文件中获取了counts和tpm表达矩阵，保存所需表达矩阵和分组信息，接着就可以用这些数据进行下游各类分析啦

参考资料

Ensembl_id转换与gene symbol基因名去重复的两种方法 - 简书 (jianshu.com)

获取基因有效长度的N种方法Counts FPKM RPKM TPM 的转化

本实战教程基于以下生信技能树分享的视频：

【生信技能树】转录组测序数据分析_哔哩哔哩_bilibili

【生信技能树】GEO数据库挖掘_哔哩哔哩_bilibili

文末友情宣传

强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：

• 数据挖掘（GEO,TCGA,单细胞）2022年5~6月场，快速了解一些生物信息学应用图表
• 生信入门课-2022年5~6月场，你的生物信息学第一课