宏基因组TPM定量

（建议学习一下，至少了解清楚TPM计算过程，否则得到比对结果后自己不知道怎么算TPM）

简单来说，TPM定量分为三步：

1. 基因长度标准化。测序reeads数除以基因长度length，得到RPK

2. 进行百万转换。因为计算的单位为M

3. 测序深度标准化。reads数除以总reads数即总RPK

宏基因组的定量用TPM是个不错的选择。

进行宏基因组的TPM定量，需要文件预测的CDS序列文件和对应的fastq（cleandata）。

软件安装：BWA和Samtools

#傻瓜式安装conda
#samtools
conda install -c bioconda samtools

#bwa 
conda install -c bioconda bwa
#如果你会自己编译环境可以不用conda安装，反正能装上并且能使用就行

以bwa构建CDS序列索引进行比对

#构建索引
bwa index Unigene.fasta

#例如我有5个样本A,B,C,D,E
for i in {A,B,C,D,E} \
do bwa mem -k 19 – t 24 Unigene.fasta \        
${i}_clean_R1.fastq ${i}_clean_R2.fastq > ./${i}.sam \ #压缩的fastq文件也可以，即后缀为gz
done

samtools提取比对到clean上的reads数

#Samtools转换sam文件为bam文件 1.3版本前 单个样本举例，多样本for循环
samtools view -bS sample.sam > sample.bam
samtools sort sample.bam > sample_sort.bam
samtools index sample_sort.bam

#1.3版本后（sort将sam转化为bam与排序同时进行）
samtools sort sample.sam > sample_sort.bam
samtools index sample_sort.bam

#获取每个ORF比对的read数
for i in {A,B,C,D,E} \
do samtools idxstats ${i}_sort.bam > ${i}_mapped.txt \
done
#注意，这一步之前需要经过sort和index

获得的txt结果有四列，从左至右分别是：基因名、基因长度、mapped_read、unmapped_read

 因为输出的文件不含表头，我们手动添加一个

for i in {A,B,C,D,E} \
do \
sed -i "1 i GeneID\t\length\tmapped_read\tunmapped_read " ${i}_mapped.txt \
done

 因为TPM计算是不需要unmapped_read的，所以只保留前面三列

for i in {A,B,C,D,E} \
do \
cut -f 1-3 ${i}_read.txt > ${i}_read_cut.txt \
done

这时候的文件就是我们需要的了，只需根据公式来计算即可以获得TPM ，这里我贴上用R语言和python计算的代码，比较简单。

R的脚本：

#读取文件
df <- read.delim("BJS_read.txt",header = T,row.names = 1)
#截取需要的部分
df_cut <- df[,1:2]
#计算RPK
df_cut$RPK <- df_cut$mapped_read*1000/df_cut$length
n <- nrow(df_cut)
result <- df_cut[-n,]#文件最后一行存在一个*行，含有NULL，需要去除掉才能计算，否则结果都为NULL
TotalRPK <- sum(result$RPK)
result$TPM <- (result$RPK*10e6)/TotalRPK
write.csv(result,file = "BJS_TPM.csv")

 python：

import pandas as pd
import argparse

# 命令行参数解析
parser = argparse.ArgumentParser()
parser.add_argument('inputfile', help='Input file name')
parser.add_argument('outputfile', help='Output file name')
args = parser.parse_args()

# 读取数据
count = pd.read_table(args.inputfile, index_col=0)

# 过滤不需要的行
count = count[count.index != '*']

# 计算 RPK 和 TPM
count['RPK'] = (count['mapped_read'] * 1000) / count['length']
total_rpk = count['RPK'].sum()
count['TPM'] = (count['mapped_read'] * 10e6) / total_rpk

# 选择需要的列
result = count[['RPK', 'TPM']]

# 保存结果
result.to_csv(args.outputfile, index=True, sep="\t")

两个的结果都一样，需要注意的是由sam文件到最后的txt文件的最后一行会含有一个*的行，代表的是未匹配到任何CDS的reads，需要在计算中将其去除。

最终会得到每个样本的TPM文件，可以使用python或者R中的merge函数将所有的样本TPM拼接到一起，得到一个最终的表格。

扫码关注微生物多组学公众号，后期会更新更多的组学干货。您的关注使我们最大的鼓励。