发布于 由Wolf_rabbit

转录组表达量分析（Counts to TPM）

此步骤计算量较小，可以用个人PC完成

有多个文件夹的结果，可以尝试用以下方法解决，单一文件夹直接下载count.txt文件即可

1、将各文件夹的txt文件统一复制到Cal_Expression文件夹内，方便后续下载到电脑。

#!/bin/bash

# 定义要循环的count.txt文件现存的路径列表

paths=(

        "/home/hztext/dai/Ana_trans/bam_file/sorted/Tc_Female"

        "/home/hztext/dai/Ana_trans/bam_file/sorted/Tc_Fruit"

        "/home/hztext/dai/Ana_trans/bam_file/sorted/Tc_Male"

        "/home/hztext/dai/Ana_trans/bam_file/sorted/Tc_Petals"

        "/home/hztext/dai/Ana_trans/bam_file/sorted/Tc_VegOr"

)

 

# 循环遍历路径列表

for path in "${paths[@]}"; do

  # 在每个路径下执行您的命令

  cd "$path"  # 切换到当前路径

 

# 遍历每个 Counts.txt 文件

for counts_file in "$path"/*.txt; do

#将$counts_file文件复制到/home/hztext/dai/Ana_trans/Cal_Expression路径中。

  cp $counts_file /home/hztext/dai/Ana_trans/Cal_Expression

done

done

2、利用WinSCP软件将所有count.txt文件下载到个人电脑，文件夹名定义为Counts

3、新建一个Rproject，将这些counts文件复制到Rproject目录下

（1）新建一个R程序，并将下列代码复制其中（无需修改，直接运行）

#清空当前环境中的所有变量

rm(list = ls())

# 定义 counts 表格所在的文件夹路径

counts_folder <- "./Counts/"

 

# 获取 counts 表格文件列表

count_file <- list.files(counts_folder, pattern = "\\.txt$", full.names = TRUE)

 

# 循环遍历每个 counts 表格进行转换

for (count_file in count_file) {

# 读取 featureCounts 结果文件

counts_data <- read.table(count_file, header = TRUE, skip = 1)

# 简化表头信息

sample_name <- gsub("./Counts/","",sub("(.*)_hisat2.counts\\.txt", "\\1", count_file))

names(counts_data)[ncol(counts_data)]<-sample_name

# 提取基因长度信息

gene_lengths <- counts_data$Length

# 计算每个样本的总reads数

sample_totals <- sum(counts_data[,ncol(counts_data)])

# 计算 TPM

tpm_data <- counts_data[, ncol(counts_data)] / (gene_lengths / 1000) / (sample_totals / 1e6)

# 添加基因名列

tpm_data <- cbind(counts_data[, 1], tpm_data)

#对tpm_data文件表头进行重命名

colnames(tpm_data)<-c("Geneid",paste0(gsub("_counts","",sample_name),"_tpm"))

head(tpm_data)

#定义输出文件名

tpm_file <- gsub("_counts","",paste0("./TPM/",sample_name,"_tpm_output.txt"))

# 将结果保存为 TPM 文件

write.table(tpm_data, file = tpm_file, sep = "\t", 

            quote = FALSE, col.names = TRUE, row.names = FALSE)

}

（2）新建合并的R程序，对输出的TPM文件进行合并（无需修改，直接运行）

# 指定要合并的文件路径

TPMs_folder <- "./TPM/"

# 创建一个空的数据框用于存储合并结果

merged_data <- data.frame()

# 获取 TPMs 表格文件列表

TPM_file <- list.files(TPMs_folder, pattern = "\\.txt$", full.names = TRUE)

# 循环读取并合并文件

for (TPM_file in TPM_file) {

  # 读取当前文件

  current_data <- read.table(TPM_file, header = TRUE, stringsAsFactors = FALSE)

  

  # 判断是否为第一个文件，如果是，则直接赋值给合并结果

  if (length(merged_data) == 0) {

    merged_data <- current_data

  } else {

    # 合并当前文件到已有结果

    merged_data <- merge(merged_data, current_data, by = "Geneid", all = TRUE)

  }

}

# 输出合并结果

write.table(merged_data, file = "all_TPM.txt", sep = "\t", 

            quote = FALSE, col.names = TRUE, row.names = FALSE)

此步骤将获得TPM值，与FPKM值有些许差别但不大，可以自行选择，只需修改计算公式即可。

 

Loading