bulk-RNA seq测序数据分析流程

假如有bulk-RNA测序的数据：TH1，TH2，TH3三个重复（实验组），TW1，TW2，TW3三个重复（对照组）

• 准备工作

需要安装的软件（如FastQC、Trimmomatic、HISAT2、StringTie、samtools）
conda install -c bioconda fastqc
conda install trimmomatic
conda install -c bioconda stringtie 
conda install -c bioconda hisat2
conda install samtools

1. 质量控制

首先，需要对原始的测序数据（FASTQ格式）进行质量控制。通常使用FastQC进行初步的质量检查，然后使用Trimmomatic或者其他工具进行质量修剪。

# 质量检查
fastqc *.fq.gz

# 质量修剪
for sample in TH1 TH2 TH3 TW1 TW2 TW3
do
    trimmomatic PE ${sample}.R1.fq.gz ${sample}.R2.fq.gz \
    ${sample}.R1.trim.fq.gz ${sample}.R1.untrim.fq.gz \
    ${sample}.R2.trim.fq.gz ${sample}.R2.untrim.fq.gz \
    SLIDINGWINDOW:4:20 MINLEN:25
done

2.读段比对

使用比对工具（如HISAT2, STAR等）将质量修剪后的读段比对到参考基因组。

# 假设参考基因组索引为genome_index
for sample in TH1 TH2 TH3 TW1 TW2 TW3
do
    hisat2 -x genome_index -1 ${sample}.R1.trim.fq.gz -2 ${sample}.R2.trim.fq.gz \
    -S ${sample}.sam
done

3.转录组组装和表达量估计

使用工具如StringTie或Cufflinks对SAM/BAM文件进行转录组组装，生成转录本。

# 转换SAM为BAM并排序
for sample in TH1 TH2 TH3 TW1 TW2 TW3
do
    samtools view -bS ${sample}.sam | samtools sort -o ${sample}.sorted.bam
done

# 使用StringTie进行转录组组装
for sample in TH1 TH2 TH3 TW1 TW2 TW3
do
    stringtie ${sample}.sorted.bam -o ${sample}.gtf
done

4.差异表达分析

最后，使用DESeq2或edgeR等工具进行差异表达分析，以识别实验组和对照组之间差异表达的基因。

library(DESeq2)

# 假设countData为一个矩阵，行为基因，列为样本，包含基因的计数
# colData为DataFrame，包含样本信息
# 这里需要根据实际情况准备countData和colData

dds <- DESeqDataSetFromMatrix(countData = countData,
                              colData = colData,
                              design = ~ condition)

dds <- DESeq(dds)
res <- results(dds)

# 保存结果
write.csv(as.data.frame(res), file="DESeq2_results.csv")