Linux的RNA-seq分析（stringtie-FPKM和HTseq计数）

本文链接：https://blog.csdn.net/weixin_43927366/article/details/136721359

一、stringtie计算基因和转录本的FPKM

1.1转录组组装

#确认安装了stringtie
#确认产生了排序的bam文件，sorted.bam
#提前安装prepDE.py，可以对stringtie输出文件进行定量
 wget http://ccb.jhu.edu/software/stringtie/dl/prepDE.py3
 将prepDE.py3放入./home/lumino/tool

#1.转录本组装
#在./SRR3418005/FPKM文件夹中输出转录本组装结果
 mkdir ./SRR3418005/FPKM
#转录本组装，这里关注新转录本，使用之前merged.gtf
stringtie -e -p 8 -G ./stringtie-merged.gtf.gtf -A ./SRR3418005/FPKM/SRR3418005.gene.gtf -o ./SRR3418005/FPKM/SRR3418005.transcripts.gtf ./SRR3418005.sorted.bam

-e：这个选项告诉 stringtie 不仅要输出组装得到的转录本，还要输出那些表达量较低、可能不被认为是完整转录本的“额外”转录本片段。这有助于捕获那些可能只部分被测序覆盖的转录本区域。
-p 8：这个选项指定了 stringtie 在执行时应使用的线程数。在这里，-p 8 表示使用两个线程，这可以加速计算过程，特别是在处理大型数据集时。
-G ./stringtie-merged.gtf：这里 stringtie-merged.gtf，是所有样本转录本组装并整合的gtf。这个文件包含了所有样本的基因和转录本的结构信息。如果不关注新转录本，可以使用基因组注释文件tair10.gtf
-A SRR3418005.gene.gtf：这个选项指定了一个输出文件的路径和名称，该文件将包含基于转录本组装结果推导出的基因级别的丰度信息。在这里，输出文件名为 SRR3418005.gene.gtf。
-o /SRR3418005/FPKM/SRR3418005.transcripts.gtf：这个选项指定了输出转录本组装结果的GTF文件的路径和名称，且包含转录本水平。输出文件名为 SRR3418005.transcripts.gtf。
../SRR3418005.sorted.bam：这是你的输入文件，一个排序后的BAM文件。它包含了RNA-Seq测序数据比对到参考基因组的结果。

1.2计算stringtie的FPKM

#2.FPKM计算
#在/SRR新建一个文件夹FPKM
#产生的SRR3418005.transcripts.gtf等文件集中到同一个文件夹
#新建一个txt
touch  sample_list.txt
#txt里按 SA01 ./SRR3418005.transcription.gtf每一行写样本名和stringtie输出的transcription.gtf
#使用python输出带所有样本的表达矩阵
python /home/lumino/tool/prepDE.py3 -i sample_list.txt -t transcripts_count.csv

python：这是启动Python解释器的命令，用于执行后面的脚本。
/home/lumino/tool/prepDE.py3：这是脚本的完整路径，它告诉Python解释器在哪里可以找到并执行prepDE.py3脚本。
-i sample_list.txt：这是传递给prepDE.py3脚本的一个参数，-i是--input的缩写，用于指定一个输入文件。在这个例子中，输入文件是sample_list.txt，这个文件应该包含了样本ID和它们对应的GTF文件路径。脚本将使用这些信息来处理每个样本的数据。
-t transcripts_count.csv：这是传递给脚本的另一个参数，-t是--transcript_count_matrix的缩写，用于指定转录本计数矩阵的输出文件名。在这个例子中，脚本将生成的转录本计数矩阵保存为transcripts_count.csv文件。

sample_list.txt按如下写

二、使用HTseq-count计算每个基因的reads数

由于同源性或序列重复，某些reads被mapping到若干个基因组位置，这是多作图现象（multi mapping），HTseq完全忽略这些多重读段。如下图，？代表reads未被计数到白色基因组。

确保同目录下有排序的bam文件以及samtolls index 过的bai文件，以及基因组注释文件gtf（ensmbl上下载的）

htseq-count -s no -r pos -f bam -i gene_id SRR3418005.BAM tair10.gtf > SRR3418005_count.txt

-s no 表示单端测序数据（如果你的数据是成对的双端测序数据，你需要根据实际情况选择yes或者提供配对信息）。
-r pos 指定读取比对到参考序列的正链。
-f bam 指定输入文件格式为BAM。
-i gene_id 是一个可选参数，它告诉htseq-count使用GTF/GFF文件中的哪个属性作为特征ID。对于某些GTF文件，可能需要使用transcript_id而不是gene_id，这取决于你的GTF文件的具体内容。
SRR3418005.BAM 是你的BAM格式的测序数据文件。
tair10.gtf 是你的参考注释文件，这里假设它是GTF或GFF格式。
> SRR3418005_count.txt 将命令的输出重定向到一个文本文件，这样你可以稍后查看或分析它。

执行这个命令后，SRR3418005_count.txt文件中将包含与tair10.tf文件中注释的每个特征匹配的reads数目。

默认采用的是并集（union）模式秒如果想改成别的模式，可以用--mode来改变。

三、合并表达矩阵

#新建一个py脚本

touch merge_count.py

笔记本上输入以下代码

import sys  
import argparse  
  
# 初始化字典和文件名列表  
mydict = {}  
filenames = []  
  
# 创建命令行参数解析器  
parser = argparse.ArgumentParser(description='Merge files and output to a specified file.')  
parser.add_argument('input_files', metavar='input_file', type=str, nargs='+',  
                     help='Input files to be merged (ending with _count.txt).')  
parser.add_argument('--output', '-o', type=str, required=True,  
                     help='Output file to write the merged data.')  
  
# 解析命令行参数  
args = parser.parse_args()  
input_files = args.input_files  
output_filename = args.output  
  
# 遍历命令行参数中的文件  
for file in input_files:  
    # 去除文件名中的_count.txt部分，并添加到文件名列表中  
    filename_without_ext = file.rsplit('_count.txt', 1)[0]  
    filenames.append(filename_without_ext)  
        
    # 读取文件的每一行并处理  
    with open(file, 'r') as f:  
        for line in f:  
            key, value = line.strip().split('\t')  
            if key in mydict:  
                mydict[key] += '\t' + value  
            else:  
                mydict[key] = value  
  
# 打开一个新文件以写入合并后的数据  
with open(output_filename, 'w') as outfile:  
    # 先写入文件名，从第二列开始（前面留一个制表符作为键的位置）  
    outfile.write('\t' + '\t'.join(filenames) + '\n')  
        
    # 遍历字典中的键值对并写入文件  
    for key, value in mydict.items():  
        outfile.write(key + '\t' + value + '\n')

最后输入python运行代码

python merge_count.py SRR3418005_count.txt SRR3418006_count.txt --output merged.count.txt

--output,指定输出文件名

四、stringtie计算的FPKM和HTSeq的差异

可以看到，Ssringtie-FPKM因为前面关注了新转录本，像AT1GO1020每个转录本被分开计算，而HTSeq只算比对的上的，总数也有减少。