使用DEseq2對RNA-seq數(shù)據(jù)進行分析，并計算FPKM和TPM。

該過程使用GenomicFeatures包獲取外顯子長度，并計算非重疊外顯子長度之和作為基因長度。選取這一因素作為基因長度是參考文章http://www.itdecent.cn/p/3c21da32d7a4

step1. 安裝并加載包：

if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager",repos = "https://mirrors.#edu.cn/CRAN/") if (!requireNamespace("GenomicFeatures", quietly = TRUE)) BiocManager::install("GenomicFeatures") if (!requireNamespace("DESeq2", quietly = TRUE)) BiocManager::install("DESeq2")

library(GenomicFeatures)

library(DESeq2)

step2. 讀取gtf文件，計算基因長度：

txdb <- makeTxDbFromGFF("84kV1_3.gtf", format = "gtf", circ_seqs = character())

? ebg <- exonsBy(txdb, by="gene")

? exon_gene_sizes <- sum(width(reduce(ebg)))

step3. 讀取DEseq2所需文件：

data_in = read.csv("mycounts1.csv", head=TRUE,row.names =1, check.names = FALSE)

countData = as.data.frame(data_in)

colData = read.csv("mymeta-hj.csv", head=TRUE)

names(colData)=c("id", "condition")

step4. 構(gòu)建DEseq2數(shù)據(jù)，并使用fpkm()計算FPKM值

dds <- DESeqDataSetFromMatrix(countData = countData,

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? colData = colData,

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? design = ~ condition)

? mcols(dds)$basepairs <- exon_gene_sizes

? fpkm_out = fpkm(dds)

? write.table(as.data.frame(fpkm_out), file="fpkm_out.txt",sep="\t",quote = FALSE)

注：這一步本來使用的rowRanges(dds)<-?ebg ，但后續(xù)分析過程中發(fā)現(xiàn)，這樣導(dǎo)入基因長度可能會出現(xiàn)counts和length不匹配的情況，即將一個基因的長度賦值給另一個基因。原因是?ebg <- exonsBy(txdb, by="gene")這一步會按照基因id排序，而countDATA是按照輸入counts文件的基因順序，rowRanges并沒有按照基因id將數(shù)據(jù)合并，而是直接合并。使用mcols(dds)$basepairs的方式可以避免這一情況，將數(shù)據(jù)按照geneid合并。

step5. FPKM轉(zhuǎn)換成TPM

fpkm<-read.table("fpkm_out.txt", sep = "\t", head=TRUE, row.names = 1, check.names = FALSE)

head(fpkm)

fpkmToTpm <- function(fpkm){

? exp(log(fpkm) - log(sum(fpkm)) + log(1e6))

}

tpms <- data.frame(apply(fpkm,2,fpkmToTpm), check.names = FALSE)

head(tpms)

colSums(tpms)

write.table(as.data.frame(tpms), file="TPM_out.txt",sep="\t",quote = FALSE)