本節(jié)概覽：

從featureCounts輸出文件中獲取counts與TPM矩陣:
讀取counts.txt構(gòu)建counts矩陣；樣品的重命名和分組；counts與TPM轉(zhuǎn)換；基因ID轉(zhuǎn)換；初步過(guò)濾低表達(dá)基因與保存counts數(shù)據(jù)

從salmon輸出文件中獲取counts與TPM矩陣:
用tximport包讀取quant.sf構(gòu)建counts與TPM矩陣；樣品的重命名和分組；初步過(guò)濾低表達(dá)基因與保存counts數(shù)據(jù)

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

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

1. 讀取counts.txt構(gòu)建counts矩陣，進(jìn)行樣品的重命名和分組

###環(huán)境設(shè)置
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")

#### 對(duì)counts進(jìn)行處理篩選得到表達(dá)矩陣 ####
a1 <- fread('./counts/counts.txt',
              header = T,data.table = F)#載入counts，第一列設(shè)置為列名
colnames(a1)
counts <- a1[,7:ncol(a1)] #截取樣本基因表達(dá)量的counts部分作為counts 
rownames(counts) <- a1$Geneid #將基因名作為行名
#更改樣品名
colnames(counts)
colnames(counts) <- gsub('/home/test/align/bam/','', #刪除樣品名前綴
                   gsub('_sorted.bam','',  colnames(counts))) #刪除樣品名后綴


#### 導(dǎo)入或構(gòu)建樣本信息,  進(jìn)行列樣品名的重命名和分組####
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)  #手動(dòng)編輯構(gòu)建樣品名和分組信息
name_list <- nlgl$name_list
colnames(counts) <- name_list #更改樣品名
group_list <- nlgl$group_list
gl <- data.frame(row.names=colnames(counts), #構(gòu)建樣品名與分組對(duì)應(yīng)的數(shù)據(jù)框
                 group_list=group_list)

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

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

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

基因表達(dá)量一般以TPM或FPKM為單位來(lái)展示，所以還需要進(jìn)行，若還想轉(zhuǎn)化為FPKM或CPM可參見(jiàn)Counts FPKM RPKM TPM 的轉(zhuǎn)化與獲取基因有效長(zhǎng)度的N種方法

#### counts，TPM轉(zhuǎn)化 ####
# 注意需要轉(zhuǎn)化的是未經(jīng)篩選的counts原始矩陣
### 從featurecounts 原始輸出文件counts.txt中提取Geneid、Length(轉(zhuǎn)錄本長(zhǎng)度)，計(jì)算tpm
geneid_efflen <- subset(a1,select = c("Geneid","Length"))
colnames(geneid_efflen) <- c("geneid","efflen")  

### 取出counts中g(shù)eneid對(duì)應(yīng)的efflen
efflen <- geneid_efflen[match(rownames(counts),
                              geneid_efflen$geneid),
                        "efflen"]

### 計(jì)算 TPM 公式
 #TPM (Transcripts Per Kilobase Million)  每千個(gè)堿基的轉(zhuǎn)錄每百萬(wàn)映射讀取的Transcripts
  counts2TPM <- function(count=count, efflength=efflen){
    RPK <- count/(efflength/1000)   #每千堿基reads (Reads Per Kilobase) 長(zhǎng)度標(biāo)準(zhǔn)化
    PMSC_rpk <- sum(RPK)/1e6        #RPK的每百萬(wàn)縮放因子 (“per million” scaling factor ) 深度標(biāo)準(zhǔn)化
    RPK/PMSC_rpk              
  }  

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

3. 基因ID轉(zhuǎn)換

若上游中采用的是UCSC的基因組和gtf注釋文件，則表達(dá)矩陣行名就是我們常見(jiàn)的gene symbol基因名；若上游采用的是gencode或ensembl基因組和gtf注釋文件，那么我們就需要將基因表達(dá)矩陣行名的Ensembl_id(gene_id)轉(zhuǎn)換為gene symbol (gene_name)了。
在轉(zhuǎn)換時(shí)經(jīng)常會(huì)出現(xiàn)多個(gè)Ensembl_id對(duì)應(yīng)與一個(gè)gene symbol的情形，此時(shí)就出現(xiàn)了重復(fù)的gene symbol。此時(shí)就需要我們?cè)谶M(jìn)行基因ID轉(zhuǎn)換去除重復(fù)的gene symbol。
下面展示轉(zhuǎn)化ID并合并所有重復(fù)symbol的方法，其他基因名去重復(fù)方法參見(jiàn)Ensembl_id轉(zhuǎn)換與gene symbol基因名去重復(fù)的兩種方法 - 簡(jiǎn)書(shū) (jianshu.com)

#合并所有重復(fù)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行名對(duì)應(yīng)的symbol
table(duplicated(symbol))  #統(tǒng)計(jì)重復(fù)基因名
  
###使用aggregate根據(jù)symbol列中的相同基因進(jìn)行合并 
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列中的相同基因進(jìn)行合并 
tpm <- column_to_rownames(tpm,'Group.1')

id轉(zhuǎn)換前

id轉(zhuǎn)換后

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

我們的數(shù)據(jù)中會(huì)有很多低表達(dá)甚至不表達(dá)的基因，在后續(xù)分析中可能會(huì)影響數(shù)據(jù)的分析判斷，因此需要對(duì)低表達(dá)的基因進(jìn)行篩除處理。篩選標(biāo)準(zhǔn)不唯一，依自己數(shù)據(jù)情況而定。在這里展示篩選出至少在重復(fù)樣本數(shù)量?jī)?nèi)的表達(dá)量counts大于1的行（基因），可以看到超過(guò)一半以上的基因都被篩掉了。

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

counts_filt <- counts[keep_feature, ] #替換counts為篩選后的基因矩陣（保留較高表達(dá)量的基因）
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)換類(lèi)型type = "salmon"、以及轉(zhuǎn)錄本與基因名gene symbol對(duì)應(yīng)關(guān)系文件（即我們之前得到的t2s_vm25_gencode.txt）就可以轉(zhuǎn)換得到各基因的定量關(guān)系了。其他步驟與操作featureCounts輸出文件類(lèi)似。
這里只展示了獲取基因表達(dá)的TPM值，如果還想了解如何獲得FPKM值請(qǐng)參考文章：獲取基因有效長(zhǎng)度的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的對(duì)應(yīng)關(guān)系文件
t2s <- fread("t2s_vm25_gencode.txt", data.table = F, header = F); head(t2s)

##找到所有quant.sf文件所在路徑  導(dǎo)入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表達(dá)矩陣
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)
 

#### 導(dǎo)入或構(gòu)建樣本信息,  進(jìn)行列重命名和分組 ####
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)


#### 初步過(guò)濾低表達(dá)基因 ####
#篩選出至少在重復(fù)樣本數(shù)量?jī)?nèi)的表達(dá)量counts大于1的行（基因）
keep_feature <- rowSums(counts > 1) >= 2               #ncol(counts)/length(table(group_list)) 
table(keep_feature)  #查看篩選情況
counts_filt <- counts[keep_feature, ] #替換counts為篩選后的基因矩陣（保留較高表達(dá)量的基因）
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')

通過(guò)以上步驟，成功從featureCounts或Salmon輸出文件中獲取了counts和tpm表達(dá)矩陣，保存所需表達(dá)矩陣和分組信息，接著就可以用這些數(shù)據(jù)進(jìn)行下游各類(lèi)分析啦

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

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九欧美,1769亚洲,黄色成人av

RNA-seq入門(mén)實(shí)戰(zhàn)（三）：從featureCounts與Salmon輸出文件獲取counts矩陣

RNA-seq入門(mén)實(shí)戰(zhàn)（三）：從featureCounts與Salmon輸出文件獲取counts矩陣

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

1. 讀取counts.txt構(gòu)建counts矩陣，進(jìn)行樣品的重命名和分組

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

3. 基因ID轉(zhuǎn)換

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

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

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九 欧美,1769亚洲,黄色成人av

RNA-seq入門(mén)實(shí)戰(zhàn)（三）：從featureCounts與Salmon輸出文件獲取counts矩陣

一、從featureCounts輸出文件中獲取counts矩陣

1. 讀取counts.txt構(gòu)建counts矩陣，進(jìn)行樣品的重命名和分組

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

3. 基因ID轉(zhuǎn)換

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

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

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九欧美,1769亚洲,黄色成人av

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

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