本节概览：
- 从featureCounts输出文件中获取counts与TPM矩阵:
读取counts.txt构建counts矩阵；样品的重命名和分组；counts与TPM转换；基因ID转换；初步过滤低表达基因与保存counts数据- 从salmon输出文件中获取counts与TPM矩阵:
用tximport包读取quant.sf构建counts与TPM矩阵；样品的重命名和分组；初步过滤低表达基因与保存counts数据
RNA-seq实战系列前期文章：
RNA-seq入门实战（一）：上游数据下载、格式转化和质控清洗
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的行（基因），可以看到超过一半以上的基因都被筛掉了。
#### 初步过滤低表达基因 ####（筛选标准不唯一、依情况而定）
#筛选出至少在重复样本数量内的表达量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表达矩阵，保存所需表达矩阵和分组信息，接着就可以用这些数据进行下游各类分析啦( • ̀ω•́ )✧
RNA-seq实战系列文章：( 2022.5.21 更新)
RNA-seq入门实战（零）：RNA-seq流程前的准备——Linux与R的环境创建
RNA-seq入门实战（一）：上游数据下载、格式转化和质控清洗
RNA-seq入门实战（二）：上游数据的比对计数——Hisat2+ featureCounts 与 Salmon
RNA-seq入门实战（三）：从featureCounts与Salmon输出文件获取counts矩阵
RNA-seq入门实战（四）：差异分析前的准备——数据检查
RNA-seq入门实战（五）：差异分析——DESeq2 edgeR limma的使用与比较
RNA-seq入门实战（六）：GO、KEGG富集分析与enrichplot超全可视化攻略
RNA-seq入门实战（七）：GSEA——基因集富集分析
RNA-seq入门实战（八）：GSVA——基因集变异分析
未完待续......
本实战教程基于以下生信技能树分享的视频
【生信技能树】转录组测序数据分析_哔哩哔哩_bilibili
【生信技能树】GEO数据库挖掘_哔哩哔哩_bilibili
网友评论