RNA-seq入门实战（三）：在R里面整理表达量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)

#### 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)

#合并所有重复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')

#### 初步过滤低表达基因 ####（筛选标准不唯一、依情况而定）
#筛选出至少在重复样本数量内的表达量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')

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')