前言

關(guān)于bulk-seq分析，你是使用fpkm？還是tpm?

這個問題，可能對于初學者來說，很是糾結(jié)。因為，你在文章中有使用fpkm，也有使用tpm？也有使用count？我們在“一文了解Count、FPKM、RPKM、TPM”推文中也詳細介紹了各參數(shù)的表示的意義。以及在關(guān)于Count，F(xiàn)PKM，TPM，RPKM等表達量的計算及轉(zhuǎn)換中也介紹了它們之間的轉(zhuǎn)換。

但是，學習的一直是在不斷向前，以及不斷糾錯的過程。我們每個人都在經(jīng)歷此過程。

bulk-seq分析中，是什么標準進行后續(xù)的分析呢？

對于這個問題，我們“社群”也有在討論過？也有很多同學在咨詢這個問題。有同學認為FPKM并不是很準確，推薦使用tpm值。但是也有同學一直在使用FPKM值，或是其它值進行篩選，如count（最原始的表達量），RPM等。

對于，最生信的同學來說，無論是count，FPKM、TPM等值，都只是一個數(shù)字，僅僅只是一個數(shù)字，“不能表示任何意義”，我們只需要整體的趨勢一樣即。因此，無論是使用fpkm，或是tpm，都可以滿足我們的需求。

但是，使用不同的表達量，最大的區(qū)別在于篩選的基因數(shù)量，這個是無法避免的。

FPKM

FPKM: FPKM的全稱為Fragments Per Kilobase Million，F(xiàn)ragments Per Kilobase of exon model per Million mapped fragments(每千個堿基的轉(zhuǎn)錄每百萬映射讀取的fragments)。通俗講，把比對到的某個基因的Fragment數(shù)目，除以基因的長度，其比值再除以所有基因的總長度。注意，這里的基因長度是指基因外顯子的總長度。

RPKM

**RPKM: **Reads Per Kilobase of exon model per Million mapped reads (每千個堿基的轉(zhuǎn)錄每百萬映射讀取的reads)；

RPM

RPM: RPM/CPM: Reads/Counts of exon model per Million mapped reads (每百萬映射讀取的reads)

公式：RPM = ExonMappedReads * 10^6 /TotalMappedReads

TPM

TPM：TPM的全稱為Transcripts per million，Transcripts Per Kilobase of exon model per Million mapped reads (每千個堿基的轉(zhuǎn)錄每百萬映射讀取的Transcripts)

轉(zhuǎn)換或計算

我們都知道，無論是fpkm，tpm，rpm值都是基于每個基因的count值與Gene Length進行計算。

計算的方法也很多，但是或出現(xiàn)：使用不同的計算方式，得到的結(jié)果并不一致。

這個問題，大家也會經(jīng)常遇到，包括我自己。到這一步，就會很糾結(jié)，也會質(zhì)疑自己，然后會“卡”很久。

木知，大家知否是這樣的？

Code

我們也基于前期分享，在分享使用Count轉(zhuǎn)換FPKM或是TPM值的代碼。

我們的代碼僅提供給大家參考。

Count To FPKM

1. R語言Code

##'@計算FPKM
data <- read.table("./01.cucumber.mRNA.gene.count.txt",
                      header = T, sep = "\t", row.names = 1)
head(data)
gene_length_kb <- data$Length /1000
expression_matrix <- data[, 2:ncol(data)]
head(expression_matrix)
# # 計算每列總Counts
total_counts <- colSums(expression_matrix)

# 計算FPKM值
fpkm_matrix <- sweep(expression_matrix, 2, total_counts / 1e6, "/") # 除以總Reads
fpkm_matrix <- sweep(fpkm_matrix, 1, gene_length_kb, "/") # 除以基因長度
fpkm_matrix[1:5,1:10]

輸入數(shù)據(jù)

2. perl腳本

use strict;
 
open A,"$ARGV[0]";
 
<A>;
my @colsum;
while(<A>){
        chomp;
        my @line=split;
        shift @line;
        shift @line;
        for(0..$#line){
                $colsum[$_]+=$line[$_];
        }
}
 
open B,"$ARGV[0]";
my $head=<B>;
chomp($head);
print "$head\n";
while(<B>){
        chomp;
        my @line=split;
        my $gene=shift @line;
        my $length=shift @line;
        print "$gene\t";
        for(my $i=0;$i<@line;$i++){
                my $fpkm=$line[$i]*1000000*1000/($colsum[$i]*$length);
                print "$fpkm\t";
        }
        print "\n";
}

運行

perl CountToFPKM.pl Input_count.txt > OutPut_FPKM.txt 

Count To TPM

##'@使用count計算tpm
count_df <- read.table("./01.cucumber.mRNA.gene.count.txt",
                       header = T, sep = "\t", row.names = 1)
head(count_df)

# TPM函數(shù)
calculate_tpm <- function(counts, lengths) {
  rpk <- counts / (lengths / 1000)  # RPK: count per kilobase
  scaling_factors <- colSums(rpk) / 1e6  # per million scaling factor
  tpm <- sweep(rpk, 2, scaling_factors, "/")  # divide each column by its scaling factor
  return(tpm)
}

# 準備數(shù)據(jù)
gene_lengths <- count_df$Length
count_matrix <- as.matrix(count_df[, -1])  # 排除長度列，僅保留 count 數(shù)據(jù)

# 計算 TPM
tpm_matrix <- calculate_tpm(count_matrix, gene_lengths)


# 查看結(jié)果
tpm_matrix[1:5,1:10]

關(guān)于Count值，直接使用featureCounts等軟件即可獲得。

OK，以上僅是自己的學習筆記，若是有不對的地方，請大家指出。

若是你有其他的方法，我們歡迎你的留言區(qū)進行留言討論！

若我們的教程對你有所幫助，請點贊+收藏+轉(zhuǎn)發(fā)，大家的支持是我們更新的動力??！

2024已離你我而去，2025加油?。?/strong>