发布于 由Wolf_rabbit

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

利用排序（sorted）后的bam文件，依次计算表达量

1、安装subread，其中包括 featureCounts 工具

conda install -c bioconda subread

利用featureCounts计算mRNA（根据基因组注释文件选择：mRNA、CDS、exon、gene等）表达量

（tips: 仅需要修改循环的路径列表（bam文件所在路径）、输出路径、基因组注释文件所在位置、需要计算的类别）

#!/bin/bash

# 定义要循环的路径列表

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"

)

#定义输出路径

output_path="/home/hztext/dai/copy_Ana_trans/counts_file"

mkdir -p "$output_path"

 

 

# 循环遍历路径列表

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

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

  folder_name=$(basename "$path")

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

#指定基因组注释文件路径和文件名

annotation_file="/home/hztext/dai/Ana_trans/Toona_cilita_var_Genome/Toona_ciliata_LG_Genome.gff"

 

# 循环遍历 BAM 文件

for bam_file in "$path"/*.bam; do

  # 提取文件名作为样本名

  sample_name=$(basename "$bam_file" .bam)

 

  # 使用 featureCounts 计算表达量。-T：指定线程数；-a：指定注释文件；-o指定输出路径（不填则为工作路径）/文件名；-p：指定双端测序；-t：指定输出的表达量是mRNA，gene或者exon等；-g：基因标识符属性（如'ID=Ts_Contig1_G00001.t1.exon1;Parent=Ts_Contig1_G00001.t1'则为ID）

  featureCounts -T 50 -a "$annotation_file" -o "$output_path"/"${sample_name}.counts.txt" -p -t mRNA -g ID "$bam_file"

done

done

此过程将获得以.counts.txt结尾的Counts文件

下一步将Counts值转化为TPM值……………

Loading