本節(jié)概覽：

1. 從featureCounts輸出文件中獲取counts與TPM矩陣:
   讀取counts.txt構建counts矩陣笛园；樣品的重命名和分組税手；counts與TPM轉(zhuǎn)換；基因ID轉(zhuǎn)換；初步過濾低表達基因與保存counts數(shù)據(jù)
2. 從salmon輸出文件中獲取counts與TPM矩陣:
   用tximport包讀取quant.sf構建counts與TPM矩陣书幕；樣品的重命名和分組匙监；初步過濾低表達基因與保存counts數(shù)據(jù)

承接上節(jié)RNA-seq入門實戰(zhàn)（二）：上游數(shù)據(jù)的比對計數(shù)——Hisat2與Salmon
之前已經(jīng)得到了featureCounts與Salmon輸出文件（counts凡橱、salmon）和基因ID轉(zhuǎn)化文件（g2s_vm25_gencode.txt、t2s_vm25_gencode.txt)亭姥。一般為了對樣品進行分組注釋我們還需要在GEO網(wǎng)站下載樣品Metadata信息表SraRunTable.txt稼钩，接下來就需要在R中對輸出結果進行操作，轉(zhuǎn)化為我們想要的基因表達counts矩陣达罗。

一坝撑、從featureCounts輸出文件中獲取counts矩陣

1. 讀取counts.txt構建counts矩陣，進行樣品的重命名和分組

###環(huán)境設置
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), #構建樣品名與分組對應的數(shù)據(jù)框
                 group_list=group_list)

這里給樣品名加上_1绍载、_2表示重復樣品，根據(jù)這兩類細胞的多能性狀態(tài)將其分組為naive和primed

fix(nlgl)編輯構建樣品名和分組信息

2. counts與TPM轉(zhuǎn)換

基因表達量一般以TPM或FPKM為單位來展示滔蝉，所以還需要進行击儡，若還想轉(zhuǎn)化為FPKM或CPM可參見Counts FPKM RPKM TPM 的轉(zhuǎn)化 與 獲取基因有效長度的N種方法

#### counts，TPM轉(zhuǎn)化 ####
# 注意需要轉(zhuǎn)化的是未經(jīng)篩選的counts原始矩陣
### 從featurecounts 原始輸出文件counts.txt中提取Geneid蝠引、Length(轉(zhuǎn)錄本長度)阳谍，計算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)  每千個堿基的轉(zhuǎn)錄每百萬映射讀取的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轉(zhuǎn)換

若上游中采用的是UCSC的基因組和gtf注釋文件，則表達矩陣行名就是我們常見的gene symbol基因名螃概；若上游采用的是gencode或ensembl基因組和gtf注釋文件矫夯，那么我們就需要將基因表達矩陣行名的Ensembl_id(gene_id)轉(zhuǎn)換為gene symbol (gene_name)了。
在轉(zhuǎn)換時經(jīng)常會出現(xiàn)多個Ensembl_id對應與一個gene symbol的情形吊洼，此時就出現(xiàn)了重復的gene symbol训貌。此時就需要我們在進行基因ID轉(zhuǎn)換去除重復的gene symbol。
下面展示轉(zhuǎn)化ID并合并所有重復symbol的方法冒窍，其他基因名去重復方法參見Ensembl_id轉(zhuǎn)換與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))  #統(tǒng)計重復基因名
  
###使用aggregate根據(jù)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轉(zhuǎn)換前

id轉(zhuǎn)換后

4. 初步過濾低表達基因與保存counts數(shù)據(jù)

我們的數(shù)據(jù)中會有很多低表達甚至不表達的基因递沪，在后續(xù)分析中可能會影響數(shù)據(jù)的分析判斷，因此需要對低表達的基因進行篩除處理综液。篩選標準不唯一款慨，依自己數(shù)據(jù)情況而定。在這里展示篩選出至少在重復樣本數(shù)量內(nèi)的表達量counts大于1的行（基因）谬莹，可以看到超過一半以上的基因都被篩掉了檩奠。

#### 初步過濾低表達基因 ####（篩選標準不唯一桩了、依情況而定）
#篩選出至少在重復樣本數(shù)量內(nèi)的表達量counts大于1的行（基因）
keep_feature <- rowSums(counts>1) >= 2
table(keep_feature)  #查看篩選情況，F(xiàn)ALSE為低表達基因數(shù)（行數(shù)）埠戳，TURE為要保留基因數(shù)
#FALSE  TRUE 
#35153 20339 

counts_filt <- counts[keep_feature, ] #替換counts為篩選后的基因矩陣（保留較高表達量的基因）
tpm_filt <- tpm[keep_feature, ]

#### 保存數(shù)據(jù) ####
counts_raw=counts #這里重新命名方便后續(xù)分析調(diào)用
counts=counts_filt
tpm=tpm_filt

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

二井誉、從salmon輸出文件中獲取counts矩陣

需要用到tximport包從salmon輸出文件中獲取counts矩陣，在tximport函數(shù)中輸入quant.sf文件路徑乞而、轉(zhuǎn)換類型type = "salmon"送悔、以及轉(zhuǎn)錄本與基因名gene symbol對應關系文件（即我們之前得到的t2s_vm25_gencode.txt）就可以轉(zhuǎn)換得到各基因的定量關系了。其他步驟與操作featureCounts輸出文件類似爪模。
這里只展示了獲取基因表達的TPM值欠啤，如果還想了解如何獲得FPKM值請參考文章：獲取基因有效長度的N種方法中第二部分內(nèi)容以及Counts FPKM RPKM TPM 的轉(zhuǎn)化。

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數(shù)取整
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)


#### 初步過濾低表達基因 ####
#篩選出至少在重復樣本數(shù)量內(nèi)的表達量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, ]


#### 保存數(shù)據(jù) ####
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表達矩陣洁段，保存所需表達矩陣和分組信息，接著就可以用這些數(shù)據(jù)進行下游各類分析啦

參考資料
Ensembl_id轉(zhuǎn)換與gene symbol基因名去重復的兩種方法 - 簡書 (jianshu.com)
獲取基因有效長度的N種方法
Counts FPKM RPKM TPM 的轉(zhuǎn)化
本實戰(zhàn)教程基于以下生信技能樹分享的視頻：
【生信技能樹】轉(zhuǎn)錄組測序數(shù)據(jù)分析_嗶哩嗶哩_bilibili
【生信技能樹】GEO數(shù)據(jù)庫挖掘_嗶哩嗶哩_bilibili

RNA-seq實戰(zhàn)系列文章：
RNA-seq入門實戰(zhàn)（零）：RNA-seq流程前的準備——Linux與R的環(huán)境創(chuàng)建
RNA-seq入門實戰(zhàn)（一）：上游數(shù)據(jù)下載共郭、格式轉(zhuǎn)化和質(zhì)控清洗
RNA-seq入門實戰(zhàn)（二）：上游數(shù)據(jù)的比對計數(shù)——Hisat2+ featureCounts 與 Salmon
RNA-seq入門實戰(zhàn)（三）：從featureCounts與Salmon輸出文件獲取counts矩陣
RNA-seq入門實戰(zhàn)（四）：差異分析前的準備——數(shù)據(jù)檢查
RNA-seq入門實戰(zhàn)（五）：差異分析——DESeq2 edgeR limma的使用與比較
RNA-seq入門實戰(zhàn)（六）：GO祠丝、KEGG富集分析與enrichplot超全可視化攻略
RNA-seq入門實戰(zhàn)（七）：GSEA——基因集富集分析
RNA-seq入門實戰(zhàn)（八）：GSVA——基因集變異分析
RNA-seq入門實戰(zhàn)（九）：PPI蛋白互作網(wǎng)絡構建（上）——STRING數(shù)據(jù)庫的使用
RNA-seq入門實戰(zhàn)（十）：PPI蛋白互作網(wǎng)絡構建（下）——Cytoscape軟件的使用
RNA-seq入門實戰(zhàn)（十一）：WGCNA加權基因共表達網(wǎng)絡分析——關聯(lián)基因模塊與表型