转录组分析

一、分析基本流程

（一）质控——fastqc

（二）去接头，过滤低质量reads——Trimmomatic

（三）去除rRNA——SortMeRNA

（前3步时间较久，最好上传cleandata，直接从以下步骤开始）

（四）比对——hisat2

（五）定量——featurecounts、cufflinks

（六）差异基因分析——DESeq

二、raw data分析

（一）质控——fastqc

1、手动从官网下载安装包上传到linux：https://link.zhihu.com/?target=https%3A//www.bioinformatics.babraham.ac.uk/projects/fastqc/

2、后续解压既可用

unzip fastqc_v0.12.1.zip
./fastqc -help

设置环境变量

vim ~/.bashrc ##打开文件，Ctrl+g定位最后一行，按‘i’插入以下路径
export PATH="yourdir/FastQC/:$PATH" 
：wq  ##按esc退出，直接输入‘：wq’保存退出

source ~/.bashrc ##更新路径
echo $PATH  ##查看路径
fastqc -help  ##查看fastqc是否可用

3、运行

nohup fastqc yourdatadir*.fq.gz -o youroutputdir & 
##*.fq.gz文件不加引号！！-o的路径可以加双引号！！
##测序报告解读：https://blog.csdn.net/qq_22253901/article/details/119638317
###https://mp.weixin.qq.com/s?__biz=MzA4NDAzODkzMA==&mid=2651273374&idx=1&sn=dfac9e82a09fbb18cd8b2c8179afcad8&chksm=841edce3b36955f5ef17ae96b03c11b9d0052cd3490f1d6b91e3600a5d2eafa77ad79d4a40ff&scene=21#wechat_redirect

（二）去接头，过滤低质量reads——Trimmomatic

1、安装

wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.38.zip
unzip Trimmomatic-0.38.zip
java -jar "yoursoftwaredir/Trimmomatic-0.38/trimmomatic-0.38.jar" –h
java -jar $trimmomatic ##看是否能用

2、命令参数

java -jar "yoursoftwaredir/Trimmomatic-0.38/trimmomatic-0.38.jar" PE -threads 12 -phred33 R1.fq.gz $R2.fq.gz $R1.paired.fq.gz $R1.unpaired.fq.gz $R2.paired.fq.gz $R2.unpaired.fq.gz ILLUMINACLIP:$adapter.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36
##Trimmomatic得引用绝对路径
##参数顺序很重要，不能更改
##-threads 线程数；
##-phred33： 设置碱基的质量格式，默认的是-phred64。
##最终运行出的结果输出文件有四个，R1.paired.fq.gz 、R1.unpaired.fq.gz、 R2.paired.fq.gz 、R2.unpaired.fq.gz ，如果在接下来要进行序列比对的话用的文件只需要使用到两个paired文件。
##通常的过滤步骤如下：
  ###ILLUMINACLIP: 过滤 reads 中的 Illumina 测序接头和引物序列，并决定是否去除反向互补的 R1/R2 中的 R2。
  ###SLIDINGWINDOW: 从 reads 的 5' 端开始，进行滑窗质量过滤，切掉碱基质量平均值低于阈值的滑窗。
  ###MAXINFO: 一个自动调整的过滤选项，在保证 reads 长度的情况下尽量降低测序错误率，最大化 reads 的使用价值。
  ###LEADING: 从 reads 的开头切除质量值低于阈值的碱基。
  ###TRAILING: 从 reads 的末尾开始切除质量值低于阈值的碱基。
  ###CROP: 从 reads 的末尾切掉部分碱基使得 reads 达到指定长度。
  ###HEADCROP: 从 reads 的开头切掉指定数量的碱基。
  ###MINLEN: 如果经过剪切后 reads 的长度低于阈值则丢弃这条 reads。
  ###AVGQUAL: 如果 reads 的平均碱基质量值低于阈值则丢弃这条 reads。
  ###TOPHRED33: 将 reads 的碱基质量值体系转为 phred-33。
  ###TOPHRED64: 将 reads 的碱基质量值体系转为 phred-64。
##ILLUMINACLIP:<fastaWithAdaptersEtc>:<seed mismatches>:<palindrome clip threshold>:<simple clip threshold>:<minAdapterLength>:<keepBothReads>
##ILLUMINACLIP:TruSeq3-SE:2:30:10 #接头和引物序列在 TruSeq3-SE 中，第一步 seed 搜索允许2个碱基错配，palindrome 比对分值阈值 30，simple clip 比对分值阈值 10，palindrome 模式允许切除的最短接头序列为 8bp（默认值），palindrome 模式去除与 R1 完全反向互补的 R2（默认去除）
##详细说明：https://www.jianshu.com/p/a8935adebaae

3、批量处理

（1）由于转录命名规律不一致，准备id_list让程序读取。id_list遵循以下原则：①一行一个；②行末无回车符‘\r’

sed -n l id_list.txt ##查看文件是否有回车符
cat filename | tr -d '\r' > newfilename ##删除文件回车符
##ABA1_4_1.fq.gz
##ABA1_4_2.fq.gz
##CK1_7_1.fq.gz
##CK1_7_2.fq.gz
##CK_1_8_1.fq.gz
##CK_1_8_2.fq.gz
cat id_list.txt  | awk -F '\t' '{print $1}' | sort | uniq > IDuniq.txt
##把相同后缀_1.fq.gz和_2.fq.gz替换或删除（用notepad或excel好一些），只保留一个id
##ABA1_4
##CK1_7
##CK_1_8
##CK_1_8

（2）进行循环

-d 输入的工作路径必须完整包含最后一个“\”。

#!/bin/bash
#Usage:sh run.sh -f id_list -d work_dir
# 默认值为空
id_list=""
work_dir=""

# 使用getopts来解析命令行参数
while getopts ":f:d:" opt; do
  case $opt in
    f)
      id_list="$OPTARG"  # 将命令行参数的值赋给id_list变量
      ;;
    d)
      work_dir="$OPTARG"  # 将命令行参数的值赋给work_dir变量
      ;;
    \?)
      echo "无效的选项: -$OPTARG" >&2
      exit 1
      ;;
    :)
      echo "选项 -$OPTARG 需要参数." >&2
      exit 1
      ;;
  esac
done

# 检查id_list和work_dir是否为空
if [ -z "$id_list" ] || [ -z "$work_dir" ]; then
  echo "请指定文件ID列表 (-f 参数) 和工作路径 (-d 参数)"
  exit 1
fi

# 遍历ID列表文件中的每个ID
while IFS= read -r id; do
    # 在这里，你可以执行希望针对每个ID执行的操作
    # 构建输入和输出文件的完整路径
    input_file1="${work_dir}${id}_1.fq.gz"
    input_file2="${work_dir}${id}_2.fq.gz"
    output_paired1="${work_dir}${id}_1.paired.fq.gz"
    output_unpaired1="${work_dir}${id}_1.unpaired.fq.gz"
    output_paired2="${work_dir}${id}_2.paired.fq.gz"
    output_unpaired2="${work_dir}${id}_2.unpaired.fq.gz"

    # 执行Trimmomatic命令，以ID为文件名
    java -jar "/home/share_data1/luchh/software/Trimmomatic-0.38/trimmomatic-0.38.jar" PE -threads 12 -phred33 \
    "$input_file1" "$input_file2" \
    "$output_paired1" "$output_unpaired1" \
    "$output_paired2" "$output_unpaired2" \
    ILLUMINACLIP:$adapter.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36
done < "$id_list"

3、nohup 输出正确的日志

##TrimmomaticPE: Started with arguments:
## -threads 12 -phred33 yourdatadir/bark_4_1.fq.gz yourdatadir/bark_4_2.fq.gz yourdatadir/bark_4_1.paired.fq.gz yourdatadir/bark_4_1.unpaired.fq.gz yourdatadir/bark_4_2.paired.fq.gz yourdatadir/bark_4_2.unpaired.fq.gz ILLUMINACLIP:.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36
##java.io.FileNotFoundException: /home/share_data1/luchh/TS_trans/.fa (No such file or directory)
##	at java.io.FileInputStream.open0(Native Method)
##	at java.io.FileInputStream.open(FileInputStream.java:195)
##	at java.io.FileInputStream.<init>(FileInputStream.java:138)
##	at org.usadellab.trimmomatic.fasta.FastaParser.parse(FastaParser.java:54)
##	at org.usadellab.trimmomatic.trim.IlluminaClippingTrimmer.loadSequences(IlluminaClippingTrimmer.java:110)
##	at org.usadellab.trimmomatic.trim.IlluminaClippingTrimmer.makeIlluminaClippingTrimmer(IlluminaClippingTrimmer.java:71)
##	at org.usadellab.trimmomatic.trim.TrimmerFactory.makeTrimmer(TrimmerFactory.java:32)
##	at org.usadellab.trimmomatic.Trimmomatic.createTrimmers(Trimmomatic.java:59)
##	at org.usadellab.trimmomatic.TrimmomaticPE.run(TrimmomaticPE.java:552)
##	at org.usadellab.trimmomatic.Trimmomatic.main(Trimmomatic.java:80)
##Input Read Pairs: 23082826 Both Surviving: 22680932 (98.26%) Forward Only Surviving: 297578 (1.29%) Reverse Only Surviving: 74977 (0.32%) Dropped: 29339 (0.13%)
##TrimmomaticPE: Completed successfully

（三）去除rRNA——SortMeRNA

1、下载rRNA数据库，

wget https://github.com/biocore/sortmerna/releases/download/v4.2.0/database.tar.gz
tar -zxvf database.tar.gz

2、安装

wget https://github.com/biocore/sortmerna/releases/download/v4.2.0/sortmerna-4.2.0-Linux.sh
bash sortmerna-4.2.0-Linux.sh --help
bash sortmerna-4.2.0-Linux.sh --skip-license

3、添加环境变量（同fastqc）

4、查看SortMeRNA是否可用

sortmerna --version
sortmerna -h

5、命令参数

sortmerna --ref FILE [--ref FILE] --reads FWD_READS [--reads REV_READS] [OPTIONS] ##必要参数就两个
##--ref是参考数据库
###一个是第1步下载的：/home/share_data1/luchh/database/share_bio/database/rRNA/smr_v4.3_default_db.fasta
###一个是SortMeRNA自带的8个数据库：/home/share_data1/luchh/software/SortMeRNA/sortmerna/data/rRNA_databases/...
####silva-euk-28s-id98.fasta；silva-euk-18s-id95.fasta；silva-bac-23s-id98.fasta；silva-bac-16s-id90.fasta；silva-arc-23s-id98.fasta；silva-arc-16s-id95.fasta；rfam-5s-database-id98.fasta；rfam-5.8s-database-id98.fasta
##--reads：需要原始读数文件（FASTA/FASTQ/FASTA.GZ/FASTQ.GZ），对有成对读数的文件使用两次。
##--fastx：将对齐的读数输出到FASTA/FASTQ文件中
##--aligned：对齐读取文件的前缀 [dir/][pfx] WORKDIR/out/aligned。对齐输出的目录和文件前缀，即每个输出文件都将进入指定的目录，并有指定的前缀。适当的扩展名（fasta|fastq|blast|sam|etc）将被自动添加。dir "和 "pfx "都是可选的。应该放到其他单独文件夹，这样就不混了
##--other：用法与--aligned一样，但必须与'--fastx'一起使用，这个是比对不上的，需要用于下一步比对定量。  
##--out2: 将双端文件单独输出。！！！如果没有用--out2参数，双端测序的输出数据没有分开输出，直接合并了，所以以下的hisat分析会出现问题！得重跑！！

6、批量处理（循环脚本）

rmRNA_dir="/home/share_data1/XXX/Cdtreatment/" ##*.paired.fq.gz文件所在目录
for file in ${rmRNA_dir}*.paired.fq.gz;  ##遍历.paired.fq.gz文件并运行sortmerna
do
	filename=$(basename "$file" | cut -d_ -f1) ##使用basename命令获取文件名（包含扩展名）的前缀，例如如果文件名为"sample_paired_1_1.fq.gz"，则basename命令会返回"sample_paired_1_1"。使用cut命令获取文件名前缀中的第一个字段，即去除了下划线和后缀的样本名。如果文件名前缀为"sample_paired_1_1"，则这一步将返回"sample",cut -d_ -f1,f2得到“sample_paired”。
	filename="${filename%.*}"
	sortmerna --ref "/home/share_data1/XXX/database/share_bio/database/rRNA/smr_v4.3_default_db.fasta" \
	--reads "${rmRNA_dir}${filename}_f1.paired.fq.gz" \
	--reads "${rmRNA_dir}${filename}_r2.paired.fq.gz" \
	--aligned "/home/share_data1/XXX/xunhauntest/${filename}_aligned" \ ##aligned上的需要出去，所以单独放在一个文件夹，这个文件夹层级比存放filtered的文件靠前
	--other "/home/share_data1/XXX/Cdtreatment/rRNA/${filename}_filtered" \ ##后续比对用的文件，单独放在rRNA下文件夹
	--fastx \
	--threads 10 \
	--out2 \
	--workdir "/home/share_data1/XXX/xunhauntest/${filename}/" ; ##这个路径和aligned的一样
done

其实也可以运用trimmomatic的批量处理脚本，稍加修改就行。 

#!/bin/bash

# 默认值为空
id_list=""
work_dir=""

# 使用getopts来解析命令行参数
while getopts ":f:d:" opt; do
  case $opt in
    f)
      id_list="$OPTARG"  # 将命令行参数的值赋给id_list变量
      ;;
    d)
      work_dir="$OPTARG"  # 将命令行参数的值赋给work_dir变量
      ;;
    \?)
      echo "无效的选项: -$OPTARG" >&2
      exit 1
      ;;
    :)
      echo "选项 -$OPTARG 需要参数." >&2
      exit 1
      ;;
  esac
done

# 检查id_list和work_dir是否为空
if [ -z "$id_list" ] || [ -z "$work_dir" ]; then
  echo "请指定文件ID列表 (-f 参数) 和工作路径 (-d 参数)"
  exit 1
fi

# 遍历ID列表文件中的每个ID
while IFS= read -r id; do
    # 执行Trimmomatic命令，以ID为文件名
    sortmerna --ref "/home/share_data1/luchh/database/share_bio/database/rRNA/smr_v4.3_default_db.fasta" \
    --reads "${work_dir}${id}_1.paired.fq.gz" \
    --reads "${work_dir}${id}_2.paired.fq.gz" \
    --aligned "/home/share_data1/luchh/ab_111_Trans/rRNAalign/${id}_aligned" \
    --other "/home/share_data1/luchh/ab_111_Trans/rRNA/${id}_filtered" \
    --fastx \
    --threads 10 \
    --out2 \
    --workdir "/home/share_data1/luchh/ab_111_Trans/rRNAalign/${id}/" ;
done < "$id_list"

正确运行nohup输出日志如下，最后一行有完整命令：

##以上处理后得到的数据与公司给的cleandata数据大小差异较大，自己衡量用哪一个数据。

三、cleandata分析

（一）比对——hisat2

1、安装教程：https://zhuanlan.zhihu.com/p/451939113
2、添加环境变量
3、索引建立：
（1）网站下载：https://link.zhihu.com/?target=http%3A//daehwankimlab.github.io/hisat2/download/

（2）自己构建

①将gff转成gtf

##下载gffread-master.zip
https://github.com/gpertea/gffread.git
##解压
unzip gffread-master.zip
##下载gclib
cd /home/share_data1/luchh/software/gffread-master/
https://github.com/gpertea/gclib.git
make
##修改"/home/share_data1/luchh/software/gffread-master/Makefile"，第一行加入gclib路径
GCLDIR := /home/share_data1/luchh/software/gffread-master/gclib/
##安装gffread
make release
##添加环境变量
vim ~/.bashrc
source ~/.bashrc
##查看是否添加成功
gffread -h
##转化gtf
gffread yourgenome.gff3 -T -o yourgenome.gtf

②用hisat2软件内置的脚本获取转录组和剪切位点

extract_exons.py yourgenome.gtf > yourgenome.exon
extract_splice_sites.py yourgenome.gtf > yourgenome.ss

③建索引

hisat2-build -f yourgenome.fa --ss yourgenome.ss --exon yourgenome.exon -p 7 ./yourgenome

4、准备IDlist进行比对循环

--dta-cufflinks：后续运用cufflinkks比对时需要添加该参数

#!/bin/bash
#Usage：sh script.sh -f id_list -d yourdatadir
#在倒数第二行代码设置output，不设置的话就输出在工作路径

# 默认值为空
id_list=""
work_dir=""

# 使用getopts来解析命令行参数
while getopts ":f:d:" opt; do
  case $opt in
    f)
      id_list="$OPTARG"  # 将命令行参数的值赋给id_list变量
      ;;
    d)
      work_dir="$OPTARG"  # 将命令行参数的值赋给work_dir变量
      ;;
    \?)
      echo "无效的选项: -$OPTARG" >&2
      exit 1
      ;;
    :)
      echo "选项 -$OPTARG 需要参数." >&2
      exit 1
      ;;
  esac
done

# 检查id_list和work_dir是否为空
if [ -z "$id_list" ] || [ -z "$work_dir" ]; then
  echo "请指定文件ID列表 (-f 参数) 和工作路径 (-d 参数)"
  exit 1
fi

# 进入工作路径
cd "$work_dir" || exit 1

# 遍历ID列表文件中的每个ID
while IFS= read -r id; do
    # 检查是否存在对应的文件
    if [ ! -f "${id}_1.clean.fq.gz" ] || [ ! -f "${id}_2.clean.fq.gz" ]; then
      echo "文件 ${id}_1.clean.fq.gz 或 ${id}_2.clean.fq.gz 不存在"
      continue  # 跳过当前ID并处理下一个
    fi

    # 执行hisat2命令，以ID为文件名
    hisat2 --dta -t -p 10 -x yourindexname --min-intronlen 20 --max-intronlen 4000 --rna-strandness RF \
    -1 "yourdatadir/${id}_1.clean.fq.gz" -2 "yourdatadir/${id}_2.clean.fq.gz" --dta-cufflinks  -S "youroutputdir/${id}.sam" 
done < "$id_list"

5、文件格式转化，sam to bam

#!/bin/bash
#Usage：sh script.sh -f id_list -d yourdatadir
#在倒数第二行代码设置output，不设置的话就输出在工作路径

# 默认值为空
id_list=""
work_dir=""

# 使用getopts来解析命令行参数
while getopts ":f:d:" opt; do
  case $opt in
    f)
      id_list="$OPTARG"  # 将命令行参数的值赋给id_list变量
      ;;
    d)
      work_dir="$OPTARG"  # 将命令行参数的值赋给work_dir变量
      ;;
    \?)
      echo "无效的选项: -$OPTARG" >&2
      exit 1
      ;;
    :)
      echo "选项 -$OPTARG 需要参数." >&2
      exit 1
      ;;
  esac
done

# 检查id_list和work_dir是否为空
if [ -z "$id_list" ] || [ -z "$work_dir" ]; then
  echo "请指定文件ID列表 (-f 参数) 和工作路径 (-d 参数)"
  exit 1
fi

# 进入工作路径
cd "$work_dir" || exit 1

# 遍历ID列表文件中的每个ID
while IFS= read -r id; do
    # 检查是否存在对应的文件
    if [ ! -f "${id}.sam" ]; then
      echo "文件 ${id}.sam 不存在"
      continue  # 跳过当前ID并处理下一个
    fi

    # 执行samtools命令，以ID为文件名 
    if ! samtools view -Su "${id}.sam" | samtools sort - -@ 4 -o "yourouputdir/${id}.sorted.bam"; then
      echo "samtools命令失败: ${id}"
      # 可以添加其他错误处理逻辑，如日志记录或退出
    fi
done < "$id_list"

（二）基因定量cufflinks，bam to cxb

1、hisat中添加--dta-cufflinks参数后，使用cufflink中cuffquant进行定量

（1）创建保存结果的文件夹

        -o：指定结果输出目录：包含结果文件abundances.cxb，结果文件名称一致，在同一个目录下会被覆盖，所以需要先用以下代码根据id_list创建文件夹来区别结果。

#!/bin/bash
#Usage：sh script.sh -f id_list.txt

# 默认ID列表为空
id_list=()

# 解析命令行参数
while getopts "f:" opt; do
  case $opt in
    f)
      # 从文件中读取ID列表
      id_file="$OPTARG"
      if [ -f "$id_file" ]; then
        id_list=($(cat "$id_file"))
      else
        echo "文件 '$id_file' 不存在."
        exit 1
      fi
      ;;
    \?)
      echo "无效选项: -$OPTARG" >&2
      exit 1
      ;;
  esac
done

# 获取当前路径
current_path=$(pwd)

# 遍历ID列表，创建文件夹
for folder_name in "${id_list[@]}"; do
    # 拼接文件夹路径
    folder_path="$current_path/$folder_name"

    # 检查文件夹是否已经存在，如果不存在则创建
    if [ ! -d "$folder_path" ]; then
        mkdir "$folder_path"
        echo "已创建文件夹: $folder_name"
    else
        echo "文件夹已存在: $folder_name"
    fi
done

（2）cuffquant运行

-p ：指定线程数

-b：参考基因组fasta文件

-u ：之前建立hisat2时转化的gtf文件

最后加上Hisat2产生的该样本的比对结果文件：${id}.sorted.bam

nohup cuffquant -o youroutputcreatedbyID/${id} -b "yourdatadir/genome.fa" -u "yourindexdir/genome.gtf" "yourouputdir/${id}.sorted.bam" &

结果文件abundances.cxb可用于cuffdiff等分析

（3）标准化，cxb to FPKM

-L :按顺序把表头添加，与后续cxb文件顺序对应

cuffnorm --no-update-check -o youroutputdir/cuffnorm/ -p 32 -L BL_3,BL_4,BL_5,BW_3,BW_4,BW_5 "yourindexdir/Pto.gtf" "yourcuffquantoutputdir/BL_3/abundances.cxb" "yourcuffquantoutputdir/BL_4/abundances.cxb" "yourcuffquantoutputdir/BL_5/abundances.cxb" "yourcuffquantoutputdir/BW_3/abundances.cxb" "yourcuffquantoutputdir/BW_4/abundances.cxb" "yourcuffquantoutputdir/BW_5/abundances.cxb"

以下输出结果主要关注genes.fpkm_table

（三）基因定量featurecounts

nohup featureCounts -p --countReadPairs -T 25 -o ***.count -a "/home/share_data1/luchh/database/share_bio/database/genome/P.trichocarpa/Ptrichocarpa_533_v4.1.gene_exons.gtf" -g gene_id "/home/share_data1/luchh/Cdtreatment/Potri_index/sort/CRR129703.sorted.bam" "/home/share_data1/luchh/Cdtreatment/Potri_index/sort/CRR129702.sorted.bam" "/home/share_data1/luchh/Cdtreatment/Potri_index/sort/CRR129701.sorted.bam" "/home/share_data1/luchh/Cdtreatment/Potri_index/sort/CRR129700.sorted.bam" "/home/share_data1/luchh/Cdtreatment/Potri_index/sort/CRR129699.sorted.bam" "/home/share_data1/luchh/Cdtreatment/Potri_index/sort/CRR129698.sorted.bam" & 
##-o 结果命名最好是.count结尾！！
##counts转化为FPKM：https://zhuanlan.zhihu.com/p/569512866
##批量处理
##for i in `ls *.bam`
##do
##i=${i%.bam}
##featureCounts -T 16 -f -t gene -g gene_id -a ../path/to/reference.gtf -o ${i}.count ${i}.bam
##done