[]

基因异常的表现形式在转录本层面通常表现为数量上的增加，这一现象通常通过以下三种途径实现：基因的复制增多、表达水平的提升，以及基因间的融合。

此定义如下，源自别处摘录，该总结尤为出色，特此引用以供展示，（参考来源：https://mp..qq.com/s/UQ）。

在生物信息学领域，转录组数据是最常被使用的，转录组的高表达背后dnastar key，按照我现在的理解，主要有三种逻辑。因此，从基因组角度探究转录组异常的原因，可以是一种很有趣的叙事方式。通过这种筛选得到的key，在实验中更容易得到证实。毕竟，对于细胞来说dnastar key，转录组和蛋白组相对容易发生变化，而真正能够遗传下去的，往往是基因组和表观组层面的。近期，我聆听了江涛院士关于其团队在融合基因研究方面的成果介绍，对其技术价值有了深刻认识。因此，我编写了一个脚本，它能够轻松地对获得的fastq.gz文件进行融合基因的鉴定与统计。只需调整相关参数，问题便能迎刃而解，整个过程仅需几行代码即可完成。

由此开始正文。

该脚本的设计核心在于，在获取到原始的f*q.gz文件后，首先执行质控流程，接着进行过滤和去除接头序列的操作，随后利用star软件包进行一体化分析，并最终完成融合基因的统计分析。

使用conda命令创建一个名为fusion的新环境，指定使用bioconda和conda-forge两个渠道。conda activate rna-seq
使用conda命令，在bioconda通道中，安装名为fastqc的软件包，同时，也要安装star-fusion软件包。
使用conda工具为该流程构建一个独立环境，并为其设定名称为fusion。# 下载参考数据（这个最好下载最新的）
# 首先cd到自己存贮 参考基因组的目录mkdir gene_refcd gene_ref下载链接为https://data.broadinstitute.org/Trinity/CTAT_RESOURCE_LIB/，其中包含了GRCh38_gencode_v44_CTAT_lib_Oct292023.plug-n-play.tar.gz文件。tar -xzvf GRCh38_*.tar.gz
# 构建STAR索引启动STAR工具，进入基因组生成模式。确保使用路径 ./GRCh38_gencode_v44_CTAT_lib_Oct292023.plug-n-play/ctat_genome_lib_build_dir/ref_genome.fa.star.idx，对基因组目录进行限制，不得进行任何修改。将基因组fasta文件指定为位于"./GRCh38_gencode_v44_CTAT_lib_Oct292023.plug-n-play/ctat_genome_lib_build_dir/ref_genome.fa"的路径。运行sjdbGTFfile命令，指定文件路径为./GRCh38_gencode_v44_CTAT_lib_Oct292023.plug-n-play/ctat_genome_lib_build_dir/ref_annot.gtf，对数据进行处理。    --runThreadN 32     
 将此脚本内容保存为一个sh文件，确保其具备执行权限，并命名为 fusion_pipeline.sh。

#!/bin/bash融合基因检测流程脚本 fusion_fastp_pipeline_with_stats.sh，采用fastp工具，并包含结果统计分析功能。
usage() {    echo "用法: $0 [选项]"    echo "选项:"    echo 请指定原始数据存放的文件夹路径（该项为必填项，若未指定，系统将默认设置为当前目录下的rawdata子文件夹）。    echo "  -o  输出根目录 (默认: ./results)"    echo "  -r  参考基因组路径 (必须)"    echo "  -t  并行线程数 (默认: 16)"    echo "  -h  显示帮助信息"    exit 1}
# 默认参数RAW_DIR="./rawdata"OUT_DIR="./results"REF_GENOME=THREADS=16
# 解析参数while getopts "i:o:r:t:h" opt; do    case $opt in        i) RAW_DIR="$OPTARG" ;;        o) OUT_DIR="$OPTARG" ;;        r) REF_GENOME="$OPTARG" ;;        t) THREADS="$OPTARG" ;;        h) usage ;;        *) echo "无效参数"; usage ;;    esacdone
# 参数验证[[ -z "$REF_GENOME" ]] && { echo "错误：必须指定参考基因组路径(-r)"; exit 1; }[[ $THREADS =~ ^[0-9]+$ ]] || { echo "错误：线程数必须为整数"; exit 1; }
# 日志记录LOG_FILE="${OUT_DIR}/pipeline.log"mkdir -p "$OUT_DIR" || { echo "无法创建输出根目录 $OUT_DIR"; exit 1; }exec > >(tee -a "$LOG_FILE") 2>&1
# 质量控制与过滤process_data() {    local CLEAN_DIR="${OUT_DIR}/cleandata"    mkdir -p "$CLEAN_DIR" || { echo "无法创建 $CLEAN_DIR 目录"; exit 1; }
    # 原始数据质控    local QC_RAW_DIR="${OUT_DIR}/qc_raw"    mkdir -p "$QC_RAW_DIR" || { echo "无法创建 $QC_RAW_DIR 目录"; exit 1; }    echo "[$(date)] 原始数据质控"    fastqc -t $THREADS -o "$QC_RAW_DIR" "$RAW_DIR"/*.f*q.gz
    # 使用fastp过滤    for R1 in "$RAW_DIR"/*_R1.f*q.gz; do        local SAMPLE=$(basename "${R1%_R1*}")        local R2="${R1/_R1/_R2}"        local JSON_REPORT="${CLEAN_DIR}/${SAMPLE}_fastp.json"        local HTML_REPORT="${CLEAN_DIR}/${SAMPLE}_fastp.html"
        echo "[$(date)] 处理样本: $SAMPLE"        fastp \            --in1 "$R1" --in2 "$R2" \            --out1 "${CLEAN_DIR}/${SAMPLE}_R1.clean.fq.gz" \            --out2 "${CLEAN_DIR}/${SAMPLE}_R2.clean.fq.gz" \检测适配器针对个人使用            --cut_right \确保质量分数达到15的合格标准不合格百分比限制设定为40%。规定明确，必须满足长度要求，具体数值为36。            --thread $THREADS \            --json "$JSON_REPORT" \            --html "$HTML_REPORT"    done
    # 过滤后质控    local QC_CLEAN_DIR="${OUT_DIR}/qc_clean"    mkdir -p "$QC_CLEAN_DIR" || { echo "无法创建 $QC_CLEAN_DIR 目录"; exit 1; }    echo "[$(date)] 过滤数据质控"    fastqc -t $THREADS -o "$QC_CLEAN_DIR" "${CLEAN_DIR}"/*.clean.fq.gz}
# 融合基因检测run_fusion() {    local CLEAN_DIR="${OUT_DIR}/cleandata"    local FUSION_OUT="${OUT_DIR}/fusion_results"    mkdir -p "$FUSION_OUT" || { echo "无法创建 $FUSION_OUT 目录"; exit 1; }
    for R1 in "$CLEAN_DIR"/*_R1.clean.fq.gz; do        local SAMPLE=$(basename "${R1%_R1*}")        local R2="${R1/_R1/_R2}"
        echo "[$(date)] 检测样本: $SAMPLE"        STAR-Fusion \            --left_fq "$R1" \            --right_fq "$R2" \            --genome_lib_dir "$REF_GENOME" \            --output_dir "${FUSION_OUT}/${SAMPLE}" \            --CPU $THREADS    done}
# 结果统计generate_report() {    local FUSION_OUT="${OUT_DIR}/fusion_results"    local REPORT_DIR="${OUT_DIR}/final_report"    mkdir -p "$REPORT_DIR" || { echo "无法创建 $REPORT_DIR 目录"; exit 1; }
    # 合并所有样本的融合结果    cat "$FUSION_OUT"对星融合.fusion_predictions.tsv文件执行命令，使用awk工具进行处理。'NR == 1 || !/^#FusionName/' > "${REPORT_DIR}/merged_fusion_results.tsv"
    # 统计融合基因的数量    total_fusions=$(wc -l < "${REPORT_DIR}/merged_fusion_results.tsv")    lettotal_fusions减去1，赋值给total_fusions。# 减去表头
    # 统计每个融合基因出现的次数    awk -F'\t' 若NR大于1，则对序列中的第一个元素执行计数加一操作，完成计数后，遍历计数字典中的每一个基因，输出基因名称及其对应的计数值。 "${REPORT_DIR}/merged_fusion_results.tsv" | sort -k2 -nr > "${REPORT_DIR}/fusion_gene_counts.txt"
    # 统计支持融合的跨接reads数和片段数的总和将文件中的数据按照分隔符进行分割，并计算每个分隔符分隔出的字段的总和，赋值给变量junction_read_sum。'\t' 当NR大于1时，将$6的值加到sum上，操作结束后，输出sum的值。 "${REPORT_DIR}/merged_fusion_results.tsv")将文件内容按照特定分隔符分割，并计算包含特定分隔符的行数总和。'\t' 当NR大于1时，对sum进行累加，累加的值为$7，操作完成后，输出sum的值。 "${REPORT_DIR}/merged_fusion_results.tsv")
    # 生成综合统计报告    {        echo "总融合基因数量: $total_fusions"        echo "支持融合的跨接reads数总和: $junction_read_sum"        echo "支持融合的片段数总和: $spanning_frag_sum"        echo "融合基因出现次数统计（前10）:"        head "${REPORT_DIR}/fusion_gene_counts.txt"    } > "${REPORT_DIR}/summary_report.txt"}
# 主流程main() {    echo "====== 开始流程 $(date) ======"    process_data    run_fusion    generate_report    echo "====== 流程完成 $(date) ======"}
main
        

参数解答：

这里我指定了四个传参变量：

echo"  -i  原始数据目录 (必须, 默认: ./rawdata)"echo"  -o  输出根目录 (默认: ./results)"echo"  -r  参考基因组路径 (必须)"echo"  -t  并行线程数 (默认: 16)"

如何使用：

创建脚本.sh,

复制下面的代码：

# conda activate fusion启动命令：执行nohup指令，以bash运行名为fusion_pipeline_test.sh的脚本，并将输出结果重定向至fusion.log文件，同时将标准错误也重定向至该文件，最后在后台执行。
#这里的四个参数是必须根据自己实际情况进行填写的# 完整参数运行/export3/zhangw/Code.sum/fusion_pipeline.sh \-i 参数用于指定路径，即/export3/zhangw/mRNA/test/rawdata。将/export3/zhangw/mRNA/test/results路径下的文件进行封锁处理。请勿在/export3/zhangw/gene_ref/GRCh38_gencode_v37_CTAT_lib/GRCh38_gencode_v44_CTAT_lib_Oct292023.plug-n-play/ctat_genome_lib_build_dir目录下进行修改。    -t 48

然后运行,提交后台

nohup执行bash命令 fusion_pipeline_test.sh，并将输出重定向至文件 fusion.log。2>&1 &

输出文件结构

results/├── cleandata/             # 处理后的clean数据├── qc_raw/                # 原始数据质控报告├── qc_clean/              # 过滤数据质控报告├── fusion_results/        # 各样本融合检测结果├── final_report/│   ├── 存储了融合后的数据# 合并后的融合基因检测结果融合基因计数文件，即 fusion_gene_counts.txt。# 融合基因出现次数统计└── 摘要报告文件.txt# 综合统计报告└── pipeline.log           # 完整运行日志

如有侵权请联系删除！