传统转录组测序数据上游处理

数据处理环境

硬件	属性
系统	Ubuntu 24.04.1 LTS
CPU	13th Gen Intel® Core™ i7-13700KF × 24
内存	128.0 GB

预先安装miniconda，具体过程详见我的另一篇文章：单细胞上游分析；示例数据采用自测数据，主要参考链接如下：

1. samtools安装过程
2. RNA-Seq：从fastq到表达矩阵
3. fastp安装教程

工作目录结构大致如下：
├── 0_filename2txt_2ends.sh
├── 1_qc.sh
├── 2_hisat2_build.sh
├── 3_hisat2_compare.sh
├── 4_sam_sorted.sh
├── clean
├── hisat2
├── raw
│ └── SRR_Acc_List.txt
├── ref
│ ├── gencode.vM35.chr_patch_hapl_scaff.basic.annotation.gff3
│ ├── GRCm39.genome.fa

├── ref
└── sorted

首先生成一下项目编号文件SRR_Acc_List.txt，和fastq文件放在一起

#!/bin/bash
workd="/home/williamhan/bioproject/sah_traditional_transcriptome_2024"
cd $workd/raw

# 指定输出文件的名称
output_file="SRR_Acc_List.txt"
# 清空原文件
> $output_file
# ls ./raw >> $output_file
for file in raw/*.fastq.gz; do
	# s/old/new/ 会将文本中的 "old" 替换为 "new",可使用正则表达式
    echo "${file%%.*}" | sed 's/raw\///'
# 管道符排序并输出
done | sort -u >> $output_file

质量控制

采用国产质控软件fastp，效果非常棒，我按github上方法安装的binary文件，有点小麻烦，其他两种手段没有试过，具体详见参考链接3。

批处理脚本如下：

workd="/home/williamhan/bioproject/sah_traditional_transcriptome_2024"
cd $workd/raw

cat SRR_Acc_List.txt | while read i ; 
do
	# echo $i.raw.R1.fq.gz
	# # 单端测序
	# fastp -I SRR******.fastq -O SRR******_clean.fastq
	# 双端
	fastp -i $i.R1.raw.fastq.gz -o $i.R1.clean.fastq.gz \
	-I $i.R2.raw.fastq.gz -O $i.R2.clean.fastq.gz
done

# 把处理后文件放到clean文件夹中
cd $workd
mkdir clean
mv raw/*.clean.fastq.gz clean

构建参考序列

hisat2直接走conda安装即可，可能的坑见我的另一篇文章: Hisat2报错“libstdc++.so.6

序列文件从GENCODE - Home page (gencodegenes.org)下载，需要下载的是fasta文件和gtf文件

# hisat2构建参考序列
workd="/home/williamhan/bioproject/sah_traditional_transcriptome_2024"
cd $workd

# 新建一个ref文件夹，从gencode下载ref和gff3后构建序列
genome=ref/GRCm39.genome.fa
genome_prefix=ref
project_code=sah_traditional_transcriptome_2024
# 参数一是fasta文件地址
# 参数二是输出文件地址+前缀
hisat2-build ${genome} ${genome_prefix}/${project_code} -p 20

序列比对

# hisat2序列比对
workd="/home/williamhan/bioproject/sah_traditional_transcriptome_2024"
cd $workd

# 定义路径
path="hisat2"
# 检查路径是否存在
if [ ! -d "$path" ]; then
    # 如果路径不存在，则创建它
    mkdir -p "$path"
fi

genome_prefix=ref
project_code=sah_traditional_transcriptome_2024
# 双端
cat raw/SRR_Acc_List.txt | while read acc_id ; 
do
	hisat2 --new-summary -x ${genome_prefix}/${project_code} --threads 20 \
	-1 clean/${acc_id}.R1.clean.fastq.gz -2 clean/${acc_id}.R2.clean.fastq.gz \
	-S hisat2/${acc_id}.sam; 
done

# # 单端（未debug，修改后用！）
# for sample in `ls ~/RNASeq-analysis/data/SRR*.clean.fq.gz \
# | perl -lpe 's/SRR*.clean.fq.gz//'`; 
# do 
# 	hisat2 --new-summary -p 2 -x $genome_prefix -U $left \
# 	-S ~/RNASeq-analysis/hisat2/${sample}.sam 2> ~/RNASeq-analysis/hisat2/${sample}.err; 
# done

比对结果压缩、排序及构建索引

samtools使用conda安装

workd="/home/williamhan/bioproject/sah_traditional_transcriptome_2024"
cd $workd

# 定义路径
path="sorted"
# 检查路径是否存在
if [ ! -d "$path" ]; then
    # 如果路径不存在，则创建它
    mkdir -p "$path"
fi

sam_prefix="hisat2"
# 读取包含样本名称的文件
cat raw/SRR_Acc_List.txt | while read file ; 
do
    echo $file
    # 将 SAM 文件转换为 BAM 格式。输出文件存储在 hisat2 目录下
    samtools view -b ${sam_prefix}/${file}.sam > ${sam_prefix}/${file}.bam \
    --threads 19 --verbosity DEBUG
    # ERROR：只显示错误信息。
    # WARNING：显示错误信息和警告信息。
    # INFO：显示错误、警告和一般信息性消息。
    # DEBUG：显示所有级别的信息，包括调试信息
    # 对生成的 BAM 文件进行排序，并将排序后的文件存储在 sorted 目录下
    samtools sort ${sam_prefix}/${file}.bam > ${path}/${file}.sorted.bam \
    --threads 19 --verbosity DEBUG
    # 为排序后的 BAM 文件生成索引
    samtools index ${path}/${file}.sorted.bam \
    --threads 19
done

表达定量和标准化

multiqc和subread均使用conda安装

subread有R版本，安装过程详见我的另一篇文章：Linux安装R

# 激活环境
# conda init
# conda activate transcriptomics

# 打开featureCounts的帮助文档
# featureCounts 

# 建立工作目录
workd="/home/williamhan/bioproject/sah_traditional_transcriptome_2024"
cd $workd

# 定义路径
path="counts"
# 检查路径是否存在
if [ ! -d "$path" ]; then
    # 如果路径不存在，则创建它
    mkdir -p "$path"
fi



# 设置featureCounts所需要的gtf文件位置（注释文件）
gtf=/home/williamhan/bioproject/gencode/gencode.vM35.chr_patch_hapl_scaff.basic.annotation.gtf.gz

for file in `ls sorted/*.sorted.bam`;do
    name=$(echo "$file" | sed "s/sorted\///")
    # 运行featureCounts对所有样本进行基因水平定量
    featureCounts -T 20 -p \
    -t exon -g gene_id \
    -a $gtf -o ${workd}/counts/${name}.id.txt \
    ${file}
done
# -T 1：线程数为1
# -p：表示数据为paired-end，双末端测序数据
# -t exon表示 feature名称为exon
# -g gene_id表示meta-feature名称为gene_id(ensembl名称)
# -a $gtf 表示输入的GTF基因组注释文件；-o all.id.txt 设置输出文件类型
# *.sort.bam是被分析的bam文件。featureCounts支持通配符*

# 使用multiqc对featureCounts统计结果进行可视化。
multiqc counts/all.id.txt.summary

# 只保留all.id.txt文件的【基因名】和【样本counts】
cat counts/all.id.txt | cut -f1,7- > counts/counts.txt

输出结果为txt文件，每一行代表一个基因的信息，包括基因ID、染色体位置、起始和结束位置、链的方向（正链或负链）、基因长度以及一个与测序数据相关的值。下面是对每一列的简要解释：

1. Geneid: 基因的唯一标识符。
2. Chr: 基因所在的染色体。
3. Start: 基因在染色体上的起始位置。
4. End: 基因在染色体上的结束位置。
5. Strand: 基因的方向，“+” 表示正链，“-” 表示负链，可能包含多个值，表示基因在染色体上的多个位置。
6. Length: 基因的长度，即从起始位置到结束位置的距离，计算tpm有用。
7. sorted/L1EID0802259-A1.sorted.bam: 这是与基因相关的测序数据文件名，表示该基因的表达情况。

Geneid	Chr	Start	End	Strand	Length	sorted/L1EID0802259-A1.sorted.bam
ENSMUSG00000102693.2	chr1	3143476	3144545	+	1070	0
ENSMUSG00000064842.3	chr1	3172239	3172348	+	110	0
ENSMUSG00000051951.6	chr1;chr1;chr1	3284705;3491925;3740775	3287191;3492124;3741721	-;-;-	3634	1079
ENSMUSG00000102851.2	chr1	3322980	3323459	+	480	0
ENSMUSG00000103377.2	chr1	3435954	3438772	-	2819	1

例如，示例中第一行数据表示一个基因（ENSMUSG00000102693.2）位于染色体1上，从位置3143476到3144545，是一个正链基因，长度为1070。

使用R处理生成tpm表格（待续）