基因组注释从0开始（5）：RNA-seq数据辅助组装

20 Jan 2025 6176字 21分
CC BY 4.0 （除特别声明或转载文章外）
如果这篇博客帮助到你，可以请我喝一杯咖啡~

1. RNA-seq数据

RNA-seq数据在基因组上的比对情况能够推测出内含子位置，根据覆盖度可以推测出外显子和非编码区的边界，这边我的数据有5个不同部位的转录本，分别为leaf、LY、NY、root、steam（叶、老叶、嫩叶、根、茎）。

1

1.1. RNA-seq数据比对与整合

1.1.1. 数据准备

这边需要把各个部位的RNA-seq数据建立索引并比对到基因组，这边使用hisat-2（可以conda安装）：

conda create -n hisat2
conda activate hisat2
conda install -c bioconda hisat2

这边hisat-2比对还需要参考基因组，也就是咱们已经T2T组装完成的基因组数据经过重复序列屏蔽后的数据，有hap1_masked.fa、hap2_masked.fa，需要把这两个基因组整合到一起，作为参考基因组ref_masked.fa：

cd 2.RNA-seq
cat ../hap1_masked.fa ../hap2_masked.fa > ref_masked.fa

1.1.2. Hisat2 比对

使用hisat-2结合参考基因组将各个不同部位RNA-seq数据比对

运行以下脚本创建 Hisat2 索引：

文件名：step1_build_hisat2_index.sh

#!/bin/bash
# build_hisat2_index.sh
# 构建 HISAT2 基因组索引
# 输入文件：参考基因组 ref_masked.fa
# 输出文件：基因组索引文件 ref_masked.fa_index.*

# 输入参考基因组
genome="ref_masked.fa"

# 构建索引
hisat2-build ${genome} ${genome}_index

echo "HISAT2 基因组索引已构建完成，输出文件以 ${genome}_index 为前缀"

将5个不同部位的样本进行比对：

文件名：step2_align_samples_hisat2.sh

#!/bin/bash
# hisat2-align.sh
genome="ref_masked.fa"
index="${genome}_index"

# leaf
leaf_1="../01.RawData/leaf/leaf_1.fq.gz"
leaf_2="../01.RawData/leaf/leaf_2.fq.gz"
hisat2 -x ${index} -1 ${leaf_1} -2 ${leaf_2} -S leaf.sam -p 20

# LY
LY_1="../01.RawData/LY/LY_1.fq.gz"
LY_2="../01.RawData/LY/LY_2.fq.gz"
hisat2 -x ${index} -1 ${LY_1} -2 ${LY_2} -S LY.sam -p 20

# NY
NY_1="../01.RawData/NY/NY_1.fq.gz"
NY_2="../01.RawData/NY/NY_2.fq.gz"
hisat2 -x ${index} -1 ${NY_1} -2 ${NY_2} -S NY.sam -p 20

# root
root_1="../01.RawData/root/root_1.fq.gz"
root_2="../01.RawData/root/root_2.fq.gz"
hisat2 -x ${index} -1 ${root_1} -2 ${root_2} -S root.sam -p 20

# stem
stem_1="../01.RawData/steam/steam_1.fq.gz"
stem_2="../01.RawData/steam/steam_2.fq.gz"
hisat2 -x ${index} -1 ${stem_1} -2 ${stem_2} -S stem.sam -p 20

也可以用以下bash脚本：

文件名：step2_align_samples_hisat2.sh

#!/bin/bash
# 使用 HISAT2 对每个样本进行比对
# 输入文件：配对的 FASTQ 文件
# 输出文件：SAM 文件（每个样本单独输出）

# 定义变量
genome="ref_masked.fa"   # 参考基因组
index="${genome}_index"  # HISAT2 索引
threads=20               # 并行线程数
local="../01.RawData"    # 原始数据文件路径的前缀

# 样本列表及对应文件路径
declare -A samples
samples["leaf"]="${local}/leaf/leaf_1.fq.gz ${local}/leaf/leaf_2.fq.gz"
samples["LY"]="${local}/LY/LY_1.fq.gz ${local}/LY/LY_2.fq.gz"
samples["NY"]="${local}/NY/NY_1.fq.gz ${local}/NY/NY_2.fq.gz"
samples["root"]="${local}/root/root_1.fq.gz ${local}/root/root_2.fq.gz"
samples["stem"]="${local}/steam/steam_1.fq.gz ${local}/steam/steam_2.fq.gz"

# 对每个样本进行比对
for sample in "${!samples[@]}"; do #! 表示获取数组的所有键名。samples[@] 是关联数组 samples 的值列表。
    # 分解配对的 FASTQ 文件路径
    read -r fq1 fq2 <<< "${samples[$sample]}"
    echo "正在对样本 ${sample} 进行比对..."
    
    # 使用 HISAT2 进行比对
    hisat2 -x ${index} -1 ${fq1} -2 ${fq2} -S ${sample}.sam -p ${threads}
    
    echo "样本 ${sample} 比对完成，生成文件 ${sample}.sam"
done

echo "所有样本比对完成！"

运行此脚本后，您将生成以下 SAM 文件：

​ • leaf.sam

​ • LY.sam

​ • NY.sam

​ • root.sam

​ • stem.sam

1.2.3. 将 SAM 转换为 BAM 格式并排序

SAM 文件较大且不适合直接用于后续分析。使用 samtools 将其转换为 BAM 格式并进行排序：

文件名：step3_convert_sort_bam.sh

#!/bin/bash
# 使用 samtools 将 SAM 转换为 BAM 格式并排序
# 输入文件：SAM 文件
# 输出文件：排序后的 BAM 文件

# 样本列表
samples=("leaf" "LY" "NY" "root" "stem")

# 转换与排序
for sample in "${samples[@]}"; do
    echo "正在处理样本 ${sample}..."
    samtools view -bS ${sample}.sam | samtools sort -@ 20 -o ${sample}_sorted.bam
    # rm ${sample}.sam  # 删除原始 SAM 文件
    echo "样本 ${sample} 转换完成，生成文件 ${sample}_sorted.bam"
done

echo "所有样本以及转换BAM格式，并且以及排好序！"

运行后，将获得以下排序后的 BAM 文件：

​ • leaf_sorted.bam

​ • LY_sorted.bam

​ • NY_sorted.bam

​ • root_sorted.bam

​ • stem_sorted.bam

1.2.4. 将各个部位的bam文件合并：

使用 samtools merge 将所有样本的 BAM 文件合并:

文件名：step4_merge_bam.sh

#!/bin/bash
# 合并所有样本的 BAM 文件，生成综合 BAM 文件

threads=20  # 使用的线程数

# 找到所有排序后的 BAM 文件
bam_files=$(ls *_sorted.bam)

# 合并 BAM 文件
echo "开始合并 BAM 文件..."
samtools merge -@ ${threads} merged.bam ${bam_files}
echo "BAM 文件合并完成，生成文件 merged.bam"

运行后，将获得merged.bam 文件。

1.2.5. 转录组拼接（Trinity）

文件名：step5_trinity_assembly.sh

#!/bin/bash
# 使用 Trinity 进行转录组拼接

threads=20  # 使用的线程数
memory="100G"  # 最大内存限制
bam_file="merged.bam"  # 合并的 BAM 文件

# 运行 Trinity 拼接
echo "开始使用 Trinity 进行转录组拼接..."
Trinity --genome_guided_bam ${bam_file} \
        --genome_guided_max_intron 10000 \
        --CPU ${threads} --max_memory ${memory} \
        --output trinity_out_dir

echo "Trinity 拼接完成，结果存储在 trinity_out_dir"

1.2.6. 提取编码区（TransDecoder）

文件名：step6_transdecoder_extract.sh

#!/bin/bash
# 使用 TransDecoder 提取编码区

assembly="trinity_out_dir/Trinity.fasta"  # 拼接的转录组文件

# 运行 TransDecoder 提取编码区
echo "开始提取编码区..."
TransDecoder.LongOrfs -t ${assembly}
TransDecoder.Predict -t ${assembly}

echo "TransDecoder 提取完成，生成文件："
echo "- 长编码区：Trinity.fasta.transdecoder_dir/"
echo "- 编码区预测：Trinity.fasta.transdecoder.*"

1.2.7. 注释与修正（PASA）

文件名：step7_pasa_annotation.sh

#!/bin/bash
# 使用 PASA 对转录组进行注释与修正

# 配置文件
config_file="alignAssembly.config"
genome="ref_masked.fa"  # 基因组文件
transcripts="Trinity.fasta.transdecoder.cds"  # 编码区文件

# 运行 PASA
echo "开始使用 PASA 进行注释与修正..."
Launch_PASA_pipeline.pl -c ${config_file} \
                        -C -R -g ${genome} \
                        -t ${transcripts} \
                        --ALIGNERS blat,gmap \
                        --CPU 20

echo "PASA 注释与修正完成，生成最终训练集！"

总结：

以上，按照顺序依次运行每个脚本：

  1. 构建基因组索引：bash step1_build_hisat2_index.sh
  2. 比对样本：bash step2_align_samples_hisat2.sh
  3. 转换和排序BAM文件：bash step3_bash convert_sort_bam.sh
  4. BAM文件合并：bash step4_merge_bam.sh
  5. 转录组拼接（Trinity)：bash step5_trinity_assembly.sh
  6. 提取编码区（TransDecoder）：bash step6_transdecoder_extract.sh
  7. 注释与修正（PASA）：bash step7_pasa_annotation.sh

依赖工具：

确保已安装以下软件：

​ • HISAT2

​ • Samtools

​ • Trinity

​ • TransDecoder

​ • PASA

2. 部分依赖软件安装

差的软件可以用which查看，例如：

which hisat2

或者看conda环境中有没有创建过相关软件环境(一般分部的软件我会单独创建环境去安装)：

conda env list

2.1. HISAT2 的perl依赖冲突解决：

HISAT2 命令出现了一个与 Perl 相关的错误：loadable library and perl binaries are mismatched。这个问题通常是由于系统中存在多个版本的 Perl，而当前环境中加载的 Perl 模块或二进制文件版本与运行时的 Perl 版本不匹配。

解决办法：

运行以下命令，确认你的系统中是否存在多个版本的 Perl,同时检查 Perl 的库路径是否一致：

perl -v
perl -e 'print "@INC\n";'

输出中可以看出问题的原因是 Perl 库路径冲突。HISAT2 环境使用了 Perl 5.34，但加载的某些 Perl 库来自 Perl 5.26，这导致了版本不匹配。

你可以自己清理一下，只保存一个，不过我推荐一个方法：

PERL5LIB 是 Perl 的环境变量，用于指定加载库路径的优先级。可以通过手动设置 PERL5LIB，确保只加载 Perl 5.34.0 的路径，可以设置在hisat2环境启动时自动加载perl路径，这样就不会影响其他环境：

# mkdir -p /data/wangxingbin/SoftWare/anaconda3/envs/hisat2/etc/conda/activate.d/
touch /data/wangxingbin/SoftWare/anaconda3/envs/hisat2/etc/conda/activate.d/set_perl5lib.sh

将下面的环境配置写入到上面启动文件set_perl5lib.sh中：

export PERL5LIB=/data/wangxingbin/SoftWare/anaconda3/envs/hisat2/lib/site_perl/5.34.0:/data/wangxingbin/SoftWare/anaconda3/envs/hisat2/lib/5.34.0/x86_64-linux-thread-multi

2.2. Trinity安装

conda create -n trinity_env python=3.8
conda activate trinity_env

#依赖软件安装
conda install -c conda-forge libdb=6.2.32
conda install -c conda-forge htslib=1.21 libdeflate=1.22

#设置仓库灵活优先级，灵活安装所需软件依赖
conda config --set channel_priority flexible
# conda config --set channel_priority strict #严格优先级
conda install -c bioconda trinity=2.15.2
Trinity --version

2.3. TransDecoder安装

conda create -n transdecoder_env python=3.8
conda activate transdecoder_env

conda install -c bioconda transdecoder=5.7.1

2.4. PASA安装

conda create -n pasa_env
conda install bioconda::pasa=2.5.3