写在开头：在这一个合集中，我将详细介绍bulk RNA-seq实际数据分析的流程，有哪些注意的地方以及一些小技巧
首先是获取测序数据：如果是自己的测序数据，直接从公司给的账号中下载各个样本的RawData（fastq格式）即可。
如果想要分析公共数据，在GEO数据库输入文章中给出的登录号，如GSE184771，这些数据通常是经过质控、数据归一化、差异表达分析等处理。提交者通常会详细描述其数据的处理步骤。Series Matrix File可能会提供关于实验设计、数据处理的详细信息；supplementary file可能会提供经过处理的表达矩阵，也就是我们需要下载到R中进行分析的文件。或者进入下面的SRA Run selector，SRA数据库一般储存的是原始的测序文件，下载SRR_Acc_List1.txt，以便在Linux中下载SRA数据，可以使用它们进行完整的上游和下游分析。
#SRAToolkit 是NCBI提供的用于处理来自SRA 数据库测序数据的一个工具包。前面的常用Linux命令合集中已经介绍如何下载了
#使用其中的prefetch命令批量下载数据
#运行nohup & 用于放入后台下载，避免关闭终端导致下载中断
mkdir -p GSE184771/fastq
cd /home/GSE184771/fastq
#读取文件 SRR_Acc_List1.txt 的内容，$()将读取到的内容作为参数传递到prefetch命令中，最后-O .指定下载文件的输出目录为当前目录
nohup prefetch -O . $(<SRR_Acc_List1.txt) &
#下载完 sra 文件，我们往往需要先将 SRA文件转化为 fastq 文件才能用于后续分析。使用前面安装的 SRAToolkit 中的fastq-dump工具
vi RNA-seq.sh
# 读取包含本地 SRA 文件路径的文件SRR_Acc_List2.txt，并逐行处理
#其中fastq-dump参数：--split-files将配对reads（paired-end reads）分别输出到不同的文件中，对于一方有而一方没有的reads直接丢弃。--split-3将配对reads输出到两个文件中，并将unpaired reads输出到第三个文件中。--split-spot将双端测序分为两份,但是都放在同一个文件中。
mkdir /home/GSE184771/1.QC
cd /home/GSE184771/fastq
cat SRR_Acc_List2.txt | while read i
do
# fastq-dump 命令处理SRA 文件，并将其转换为压缩的 FASTQ 文件，输出到当前目录，fastq-dump 命令处理本地的 SRA 文件，并将其转换为压缩的 FASTQ 文件，输出到当前目录
fastq-dump --gzip --split-files $i -O /home/GSE184771/1.QC && echo "** ${i} to fastq done **"
done
bash RNA-seq.sh
fastq-dump --split-files ./SRR16060501/SRR16060501.sra -O /home/GSE184771/1.QC #或者一个一个运行不报错
现在均得到了样本的原始测序fastq文件，接下来进行质控。
#如果样本数不多，其实可以一个一个运行
#列出所有（_R1.fq.gz）的路径，并将这些路径保存到文件1中。>符号被用作输出重定向操作符，作用是将命令的输出结果写入到指定的文件中。
ls /home/RNA-seq/fastq/*_R1.fq.gz > 1
ls /home/RNA-seq/fastq/*_R2.fq.gz > 2
#使用cut命令根据/分隔符提取第5个字段（第一个字段为空，完整文件路径在第5个位置），再次使用cut根据_分隔符提取第1个字段（样本名），并将结果保存到文件0中。
ls /home/RNA-seq/fastq/*_R2.fq.gz |cut -d"/" -f 5|cut -d"_" -f 1 > 0
#paste命令合并文件0、1和2的内容（分别是样本名、_R1.fq.gz和_R2.fq.gz的路径），生成文件config
paste 0 1 2 > config
#类似于SRR16060501 /home/RNA-seq/fastq/SRR16060501_R1.fq.gz /home/RNA-seq/fastq/SRR16060501_R2.fq.gz
vi RNA-seq.sh
for i in $(cat config)
do
echo $i
arr=($i)
fq1=${arr[1]}
fq2=${arr[2]}
sample=${arr[0]}
fastqc fq1 fq2 -o /home/GSE184771/1.QC
done
cd /home/GSE184771/1.QC
#整合质控结果
multiqc ./
得到了质控结果，我以一个样本的R1为例。QC结果包括以下几部分。
其中比较重要的包括：per base sequence quality每个碱基的质量。蓝色的线表示每个碱基位置上的平均质量分数（mean quality score），一般还会有黄色的框（箱线图），表示的是每个碱基位置的质量分数分布的中位数和四分位范围（interquartile range, IQR），一般要求蓝线要在绿色区域内，箱线图在黄色区域以上。碱基质量分数与碱基出错的概率关系为Q=−10log10(P)，某个位置的碱基质量为30就是它出错的概率为0.1%
Per base sequence content在前几个碱基可以不稳定，但在后面一般是AT和CG是两条平行的线。
Sequence duplication levels峰一般在最左侧。图中显示有50%左右的reads有大于10个重复。测序数据中存在大量的重复序列，有可是PCR扩增偏好；测序深度不足，可能会导致样本中的某些序列被过度重复(当测序深度不足时，测序过程中的随机抽样误差会更显著。测序仪在读取DNA片段时是随机的，如果测序深度低，每个位置被读取的次数有限，这会导致某些片段被重复读取，而其他片段可能完全没有被读取。)；生物学因素，某些生物样本中天然存在高重复率的序列，特别是在低复杂度区域或高表达的基因区域；技术问题如接头污染或样本准备中的错误，也可能导致重复序列的增加。
Overrepresented sequences显示这条序列大量出现，占全部reads的13%左右，这很可能就是接头，或者PCR偏好性扩增造成的，一般后续需要过滤掉
Adapter content，如果有接头残留，后续一定要过滤掉
接下来就需要针对质控情况进行选择性过滤了，使用Trim-galore对数据进行过滤和质控，它将fastqc与cutadapt封装在一起。
#这个可以针对不同文件质控情况，选择不同的过滤参数进行过滤
mkdir /home/GSE184771/2.trim
cd /home/GSE184771/fastq
#使用--paired 参数时，--adapter 参数会应用于R1，--adapter2 参数会应用于R2，从而去除自定义的序列，比如overrepresented sequences。--illumina 参数会同时处理R1和R2文件中的Illumina接头序列。--quality设置质量分数阈值，低于该阈值的碱基将被修剪。--stringency 4设置修剪的严格度，至少有 4 个匹配adapter的碱基序列才能进行修剪。-j 8使用的线程数为 8。-o 输出目录。--fastqc会对修剪后的文件运行 FastQC 分析
nohup trim_galore --paired --illumina --adapter2 AAGCAGTGGTATCAACGCAGAGTACTTTTTTTTTTTTTTTTTTTTTTTTT --fastqc C1_R1.fq.gz C2_R2.fq.gz --quality 25 --stringency 4 -j 8 --fastqc -o /home/GSE184771/2.trim &
#其他根据情况进行选择性过滤
#生成的_R1_val_1.fq.gz和_R2_val_2.fq.gz: 表示经过修剪并验证配对关系后的R1和R2文件,是用于后续比对和分析的文件
Bowtie2 和 HISAT2 是用于高通量测序数据的序列比对工具。Bowtie2对短读长数据（50-100bp）比对非常高效，支持局部（local）和全局（global）比对模式，可以用于ChIP-seq、DNA-seq以及小 RNA-seq。HISAT2 专门用于处理 RNA-seq 数据；适合长读长数据，优化了比对长读长（例如 100-300bp）的能力，同时也可以高效处理基因组中的剪接事件。所以我们这里使用Hisat2进行比对
#需要准备参考物种的参考基因组和注释文件，有很多网站上的数据可供下载，如UCSC、GENCODE、Ensembl、NCBI、iGenomes等。对于大多数基因组比对和下游分析，primary assembly的 FASTA和注释 文件通常已经足够
#首先构建索引，输入参考基因组序列，构建的索引的基础名称为index_m39.genome。时间长建议后台运行
mkdir /home/GSE184771/reference
mkdir /home/GSE184771/3.align
cd /home/GSE184771/reference
#在构建 HISAT2 索引时，可以使用 --ss 和 --exon 选项来包含剪接位点和外显子信息
#使用hisat2_extract_splice_sites.py 和 hisat2_extract_exons.py 脚本从 GTF 文件中提取剪接位点和外显子信息。
hisat2_extract_splice_sites.py gencode.GRCm39_primary.annotation.gtf > m39.ss
hisat2_extract_exons.py gencode.GRCm39_primary.annotation.gtf > m39.exon
nohup hisat2-build -p 4 --ss m39.ss --exon m39.exon GRCm39.primary_genecode.genome.fa index_m39.genome &
#-x为构建的索引。-1和-2为过滤后的文件，-S为输出路径，--dta启用下游转录组分析，适用于转录组数据。--new-summary生成详细的比对统计信息
nohup hisat2 -x /home/GSE184771/reference/index_m39.genome/index_m39.genome -1 C1_R1_val_1.fq.gz -2 C1_R2_val_2.fq.gz -S /home/GSE184771/3.align/C1_aligned.sam -p 8 --dta --new-summary &
#将比对结果压缩成bam文件并根据染色体排序，以节省存储空间、提高处理效率。排序后的 BAM 文件可以被索引，索引后的 BAM 文件可以快速随机访问特定区域的数据，这对于下游的分析和可视化非常重要，很多下游分析工具要求输入的 BAM 文件是按坐标排序的（如stringtie)
cd /home/GSE184771/3.align
nohup samtools sort -@ 8 C1_aligned.sam -o C1_sorted.bam &
建立BAM索引文件
nohup samtools index ./C1_sorted.bam &
Stringtie定量，gtf文件不能是压缩文件
mkdir /home/GSE184771/4.stringtie
cd /home/GSE184771
#-G指定GTF文件，-e 参数表示只对提供的注释文件中的转录本进行表达量估算，而不进行新的转录本组装。C1_genes.list 文件是一个基因级别的表达量表格，列出了每个基因的基本信息和表达量信息。C1_count.gtf 文件是一个转录本级别的 GTF 文件，包含了转录本和外显子的信息及其表达量。
nohup stringtie ./3.align/C1_sorted.bam -G ./reference/gencode.GRCm39.annotation.gtf -e -p 8 -o ./4.stringtie/C1_count.gtf -A ./stringtie/C1_genes.list &
#由于后续R下游分析使用到的基因表达矩阵都是原始的count矩阵，不能是经过FPKM或TPM归一化，以及其他转换后的数据。所以需要从中得到原始基因表达矩阵
#汇总stringtie表达量结果，例如
vi sample_list.txt
C1 /home/zhangyina/zbw_RNA-seq/stringtie/C1_count.gtf
C2 /home/zhangyina/zbw_RNA-seq/stringtie/C2_count.gtf
C3 /home/zhangyina/zbw_RNA-seq/stringtie/C3_count.gtf
E1 /home/zhangyina/zbw_RNA-seq/stringtie/E1_count.gtf
E2 /home/zhangyina/zbw_RNA-seq/stringtie/E2_count.gtf
E3 /home/zhangyina/zbw_RNA-seq/stringtie/E3_count.gtf
#将如上内容写入sample_lsit.txt,该文件为 \t（Tab） 分隔的两列，第一列为样本名称，第二列为定量的 GTF 文件的路径，注意最后不要有空行
cat sample_lsit.txt
#从https://github.com/gpertea/stringtie，下载prepDE.py3脚本
python prepDE.py3 -i sample_list.txt
#执行完上述步骤后得到的gene_count_matrix.csv即为基因表达矩阵，可以导入到R语言用edgeR / DESeq2包进行差异表达基因分析。基因表达矩阵包含了每个基因在不同样本中的原始计数（即每个基因在每个样本中检测到的读数数量）。行表示基因，列表示样本。
#生成的transcript_count_matrix.csv 文件包含了每个转录本在不同样本中的原始计数（即每个转录本在每个样本中检测到的读数数量）。行表示转录本，列表示样本。
