欢迎批评指正

一、上游处理流程

上游处理步骤包括质量检测、质量控制、比对、定量^[2]，每一步处理数据的目的都是不同，也有相关的软件与之对应。通过质量检测，原始数据的各种问题将会呈现出来，接下来的质量控制就是为了解决原始数据的质量问题。比对是将reads比对到染色体或者基因，并生成sam或者bam文件记录比对情况以及质量。定量既是统计比对到同一位置的reads，当然并不是比对上就计数，还有一些其他的筛选条件，比如比对质量过低或者比对到多个位置的reads就不能用来计数。而免比对的定量软件kallisto和salmon不会生成sam或bam文件而直接进行定量，仅输出定量文件。以上软件都可以使用Conda安装^[3]。

二、处理工具

（一）质量检测软件

用法：fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam]

例如：fastqc 01raw_data/sample1_1_R1.fastq.gz或者fastqc 01raw_data/*

[-f fastq|bam|sam ] 指需要质量检测的文件，可以是fastq文件，也可以是.gz、bam或者sam文件。[-o output dir] 指定质量报告生成的目录，如果不指定，将会在输入文件所在目录生成质量报告html。

这份报告包括Basic Statistics、Per base sequence quality、Per tile sequence quality、Per sequence quality scores、Per base sequence content、Per sequence GC content、Per base N content、Sequence Length Distribution、Sequence Duplication Levels、Overrepresented sequences、Adapter Content。Basic Statistics给出基本的统计信息，如测序平台、序列长度、GC含量、序列总数。根据Per base sequence quality得到序列每个碱基的测序质量，在质量评分在28以上为良好，意味着该位置的碱基错误率小于1.59％。Adapter Content则给出接头的含量，如果有接头存在，则可以根据Overrepresented sequences给出的接头信息在质控阶段选择对应的接头文件^[4]。

质检报告数量不多还好，我们只需一个一个查看既可，但是如果有几十份呢？此时multiqc就派上了用场。multiqc本身没有质检的能力，只是将多个fastqc生成的html文件融合为一个html文件，这样原本要查看几十个html文件，现在只需查看一个既可，极大提高了效率。

（二）质量控制软件

质量控制阶段使用的软件是Trimmomatic。Trimmomatic 支持多线程，可以快速的去除 Reads中的接头，并根据碱基质量值修剪Reads。软件有SE和PE两种模式，分别针对单端测序和双端测序数据^[5]。

在PE模式下，输入文件为样本的Reads1和Reads2文件sample_R1.fastq.gz和sample_R2.fastq.gz，输出文件一共有四个文件，分别是样本Reads1的配对序列和未配对序列以及Reads2的配对序列和未配对序列，所以对文件命名时也需要按照输出文件的顺序，例如sample_paired_R1.clean.fastq.gz sample_unpaired_R1.clean.fastq.gz sample_paired_R2.clean.fastq.gz sample_unpaired_R2.clean.fastq。

Trimmomatic会按照所给参数的顺序处理输入文件。一般以下参数顺序进行处理。

参数ILLUMINACLIP过滤 reads 中的接头，并决定是否去除反向互补的reads；参数SLIDINGWINDOW从 reads 的 5' 端开始，进行滑窗质量过滤，切掉碱基质量平均值低于阈值的滑窗；参数LEADING从 reads 的开头切除质量值低于阈值的碱基；参数TRAILING从 reads 的末尾开始切除质量值低于阈值的碱基；参数MINLEN如果经过剪切后 reads 的长度低于阈值则丢弃这条 reads；参数AVGQUAL如果 reads 的平均碱基质量值低于阈值则丢弃这条 reads；参数TOPHRED33 将 reads 的碱基质量值体系转为 phred-33；参数TOPHRED64将 reads 的碱基质量值体系转为 phred-64，现在基本上都是phred-33。

处理结束后，得到的就是sample_paired_R1.clean.fastq.gz sample_unpaired_R1.clean.fastq.gz，sample_paired_R2.clean.fastq.gz sample_unpaired_R2.clean.fastq四个文件。同时Trimmomatic也会给出一个处理报告，显示配对reads和未配对reads所占百分比。如果显示配对reads超过90％，那么在后续步骤中只需使用sample1_paired_R1.clean.fastq.gz和sample_paired_R1.clean.fastq.gz两个文件。

（三）比对软件

主要介绍HISAT2和STAR。Tophat2团队不继续更新Tophat2而开发了HISAT2，并推荐使用HISAT2，因为其速度更快，内存占用率更小，准确率更高。而STAR更是ENCODE官方推荐的RNA-seq比对工具。无论是HISAT2还是STAR，对于Tophat2来说都有很大的优势。而且综合来讲，STAR的综合表现最好^[1]。

HISAT2和STAR的使用步骤都是先构建索引，再进行比对。用基于一定算法构建的索引文件进行比对，可以明显减少比对所需的内存和计算量，同时显著提高比对的速度以及准确率。

1、HISAT2使用

（1）构建索引

用法：hisat2-build [options]* <reference_in> <ht2_index_base>

<reference_in>是参考基因组或者参考转录本的路径，<ht2_index_base>是输出的索引文件所在的目录以及前缀。[options]*指一些其他可选参数，为-p指定线程数可以加快构建索引的速度，其他参数默认既可。

例如：hisat2-build -p 4 00ref/TAIR10.fasta 03hisat2_index/TAIR10

（2）比对

用法：hisat2 [options]* -x <ht2-idx> {-1 <m1> -2 <m2> | -U <r>} [-S <sam>]

<ht2-idx>指索引文件所在目录和前缀。<m1>和<m2>指双端测序中Reads1和Reads2文件。

<r>指单端测序文件。根据自己的测序模式选择<m1>和<m2>或者 <r>。hisat2的输出文件较单一，仅有一个sam文件，<sam>指定输出的sam文件。[options]*提供了大量可选参数，大多数选择默认既可。这里需要注意我们的reads是否具有链特异性，如果有链特异性需要对参数 --rna-strandness 进行修改。同样，可以为-p指定线程数，以提高比对速率。

示例：hisat2 --rna-strandness FR -p 4 -x 03hisat2_index/TAIR10 -1 ./02clean_data/sample1_R1_paired_clean.fq.gz -2 ./02clean_data/sample1_R2_paired_clean.fq.gz
-S 04hisat2_out/sample1.sam

使用示例比对后，仅仅会得到一个无序的sam文件，为了方便定量，还需要使用samtools软件将无序的sam文件转为有序的压缩文件bam。

示例：samtools view -b -S 04hisat2_out/sample1.sam > 05samtools_out/sample1.bam; samtools sort 05samtools_out/sample1.bam > 05samtools_out/sample1_sorted.bam

-b和-S指输入文件是sam格式，输出是bam，并用 > 将输出内容重定向到一个bam文件。samtools sort默认是按照position排序，如有需要可以加-n参数，既可按照Reads的名字排序。

2、STAR使用

（1）构建索引

示例：  STAR --runThreadN 8 \

--runMode genomeGenerate \

--genomeDir 03star_index/ \

--genomeSAindexNbases 12 \

--genomeFastaFiles 00ref/TAIR10.fasta \

--sjdbGTFfile 00ref/TAIR10.gtf \

--sjdbOverhang 149

--runThreadN指定所用的线程数为8；--runMode指定STAR要完成的动作，默认是alignReads，所以这里需要指定为genomeGenerate以生成索引文件；--genomeDir指定索引文件的生成目录；--genomeSAindexNbases指构建的索引长度，默认14，建议取10-15。该值越大会消耗越多的内存，但是检索的更快。但是对于小基因组来说，不能太大，如果索引太长就会造成索引总数少的问题，所以这里选择12。也可以通过(log2(GenomeLength)/2 - 1)计算得到；--genomeFastaFiles指定参考基因组；--sjdbGTFfile指定参考基因组的注释文件，用于构建可变剪接数据库；--sjdbOverhang指定剪接点两端的长度，默认100，建议取值 (mate_length - 1)。mate_length在FASTQC给出的报告中可以查到。

（2）比对

示例：  STAR --runThreadN 8 \

--genomeDir 03star_index/ \

--readFilesCommand zcat \

--readFilesIn ./02clean_data/sample1_R1_paired_clean.fq.gz  \

./02clean_data/sample1_R1_R2_paired_clean.fq.gz \

--outFileNamePrefix ./04star_out/sample1_R1 \

--outSAMtype BAM SortedByCoordinate \

--outBAMsortingThreadN 8 \

--quantMode TranscriptomeSAM GeneCounts

--runThreadN 指定所用线程为8；--genomeDir指索引文件所在位置；--readFilesCommand指对读入文件进行的处理，这里选择zcat是指对读入文件进行解压。--readFilesIn指定输入文件，因为我们的输入文件是.gz结尾，所以需要zcat，如果这里是没有.gz结尾，也就不需要--readFilesCommand参数了；--outFileNamePrefix指输出文件的前缀；--outSAMtype指输出文件的类型。不使用该参数则表示不输出以染色体和位置定位Reads的sam文件，使用时可以选择SAM或者BAM(SAM的压缩格式)，并且SortedByCoordinate告诉STAR此bam按照coordinate也就是position进行排序。加上这个参数，输出文件就是排过序的bam文件；--outBAMsortingThreadN 指定bam文件排序时所用的线程；--quantMode告诉STAR在定量时所采用的模式，STAR会输出所需的文件，TranscriptomeSAM 表示输出比对到转录本的sam文件；GeneCounts输出一个记录比对到各个基因上reads数的文件。

使用示例命令，将会生成7个文件。其中sample1_1Aligned.sortedByCoord.out.bam 是以bam格式记录reads比对到基因组的文件。sample1_1Log.final.out 记录了比对的统计结果。sample1_1ReadsPerGene.out.tab记录了比对到每个基因上的reads数。sample1_1Aligned.toTranscriptome.out.bam是以bam格式记录reads比对到转录本的文件。生成sample1_1Aligned.toTranscriptome.out.bam并不是必须的，比如RSEM定量时需要sample1_1Aligned.toTranscriptome.out.bam，但是featureCounts就不需要。

3、STAR和HISAT2比较

STAR的参数比HISAT2多，也就意味着STAR更加灵活，用户可以根据自己的需求灵活的改变参数，而且用户不用考虑让人头疼的链特异性问题，因为STAR可以自动判断是否为链特异性测序。而且STAR众多的输出文件可以满足不同的需求。HISAT2因为索引的优势，可以相对轻松比对跨区域的reads（可变剪切)，而Tophat2耗时久，STAR耗内存，HISAT2克服了两个的缺点^[2]。

（四）定量软件

介绍的定量软件有RSEM、featureCounts 和HTSeq-count。RSEM像是一个集成的软件，对新手甚是友好，不仅提供了定量的命令，甚至连构建表达矩阵都有相应的命令。featureCounts以快和灵活著称，半分钟既可处理2000万的reads，用户更是可以选择feature和 meta-feature进行定量。HTSeq-count同样可以选择feature和 meta-feature进行定量。

1、RSEM使用

（1）构建索引

用法：rsem-prepare-reference [options] reference_fasta_file(s) reference_name

--gtf指基因组的注释文件，并且rsem将会使用该注释文件从 reference_fasta_file(s)中提取出转录本；reference_fasta_file(s)为基因组文件；reference_name指定索引的目录以及前缀。

示例：  rsem-prepare-reference --gtf 00ref/TAIR10.gtf \

00ref/TAIR10.fasta \

05rsem_index/TAIR10

使用示例命令将会生成7个文件。

（2）定量

用法： rsem-calculate-expression [options] --alignments [--paired-end] input reference_name sample_name

Input为SAM/BAM/CRAM，也就是在比对后得到的文件；--paired-end指input是双端测序；reference_name指索引文件所在目录和前缀；sample_name为输出文件的前缀；--alignments指输入文件为SAM/BAM/CRAM。

可选择的[options]有很多，这里需要注意链特异性参数--forward-prob，如果无链特异性，不加该参数既可，如果reads1是正义链则参数为1，反之为0；--no-bam-output指不输出比对的BAM文件。

示例：rsem-calculate-expression \

--forward-prob 0 \

--paired-end \

--no-bam-output \

--alignments -p 16 -q ./04star_out/sample1_1Aligned.toTranscriptome.out.bam \

./05rsem_index/TAIR10 \

./06rsem_out/sample1_1

使用示例命令将会生成3个文件。sample1_1.genes.results是以gene_id为meta-feature的定量结果，sample1_1.isoforms.results是以transcript_id为meta-feature的定量结果。sample1_1.stat为统计文件的目录。

2、featureCounts使用

（1）定量

用法：featureCounts [options] -a <annotation_file> -o <output_file> input_file1 [input_file2]

<annotation_file>为注释文件；<output_file>指定输出文件；input_file1 [input_file2]为一系列输入文件。

[options]中需要注意的有-p指定量时将以片段计数而不是reads，此参数针对双端测序；-s同样也指有无链特异性，0为默认参数表示无特异性，1表示Reads1是正义链，2表示Reads1是反义链；-T指定线程数；-t表示计数的feature；-g表示meta-feature。

示例：featureCounts -p -s 2 -T 6 -a 00ref/TAIR10.gtf \

-o 06featurecounts_quant/sample1_1_counts.txt \

-t exon -g transcript_id \

05samtools_out/sample1_1_sorted.bam

使用示例命令将输出2个文件，sample1_1_counts.txt记录了定量结果，sample1_1_counts.txt.summary统计了定量情况。

3、HTSeq-count使用

（1）定量

用法：htseq-count [options] alignment_file gff_file

alignment_file为比对得到的文件，通常为bam；gff_file为注释文件。

[options]中需要注意的有-f 指定alignment_file的格式，为bam或sam；-r指alignment_file按什么排序，有pos和name，需要根据bam文件设置；-s指链特异性，no指无特异性，yes指reads1是正义链， reverse表示reads1是反义链；-m指定量的模式，一般选择union就可以；-a表示若比对质量小于a的值，则忽略此reads；-t指定feature；-i指定meta-feature。

示例：htseq-count -f bam -r pos -s reverse -a 10 -t exon -i transcript_id -m union ./05samtools_out/sample1_1_sorted.bam 00ref/TAIR10.gtf > ./06htseq_quant/sample1_1_counts.txt

使用示例命令可以得到1个文件。sample1_1_counts.txt记录了meta-feature的counts数，文件的末尾为定量的统计信息。

4、RSEM、featureCounts 和HTSeq-count比较

使用RSEM定量时，需要先构建索引文件，而featureCounts 和HTSeq-count用比对结果直接定量，显得方便很多，而且对于不会写提取counts脚本的用户来说，RSEM构建表达矩阵的命令同样让人惊喜。RSEM定量后的结果更加多样，有gene_id和transcript_id两类。而且count、TPM、FPKM都有，为后续差异分析提供便利。但featureCounts 和HTSeq-count只能定量所指定的meta_feature，且结果单一。featureCounts的定量速度是显而易见的快^[6]。featureCounts与HTSeq-count对待多重比对reads的态度有所不同。HTSeq-count采用全部丢弃的策略，而featureCounts更加灵活，可以通过参数-m进行处理。

（五）免比对的定量软件

主要介绍kallisto和salmon。kallisto和salmon可以免去比对步骤，直接进行定量，甚至在PC端就可以处理RNA-seq数据。

1、kallisto使用

（1）构建索引

用法：kallisto index [arguments] FASTA-files

FASTA-files为输入的转录本；[arguments]必要有-i，指定生成的索引文件；-k指定 k-mer 的长度，默认为31，必要时可以修改。

示例：kallisto index -i 03kallisto_index/TAIR10 00ref/TAIR10.transcripts.fa

使用示例命令得到一个输出文件TAIR10

（2）定量

用法：kallisto quant [arguments] FASTQ-files

FASTQ-file为样本的reads1和reads2文件。[arguments]中，-i指定索引文件；-t指定线程数；-o 指定输出文件夹；-g 指定注释文件；--rf-stranded指reads1为反义链，--fr-stranded指reads1为正义链。

示例：kallisto quant --rf-stranded -t 4 -i ./03kallisto_index/TAIR10 -o ./04kallisto_quant/sample1_1/ -g 00ref/TAIR10.gtf  02clean_data/sample1_1_R1_paired_clean.fq.gz 02clean_data/sample1_1_R2_paired_clean.fq.gz

使用示例命令可以得到三个文件，Abundance.tsv记录了定量情况

2、salmon使用

(1) 构建索引

用法：salmon index [options]

-t指转录本文件；-i指输出的索引所在目录；-p指定线程

示例：salmon index -p 12 -t 00ref/Arabidopsis_thaliana.TAIR10.cds.all.fa \

-i 03salmon_index/TAIR10

使用示例命令可以得到15个文件

(2) 定量

用法：salmon index [options]

-i 表示索引文件所在目录；-l这里需要指定测序模式，此时指定为A表示salmon自动判断；-g可以理解为meta-feature的对应关系，加上此参数，salmon还可以对gene_id进行定量，而不只是transcript_id；-1，-2分别指输入的reads1和reads2；-p指线程数；-o指输出目录。

示例：salmon quant -i 03salmon_index/TAIR10 -l A -g 1.txt \

-1 02clean_data/sample1_1_R1_paired_clean.fq.gz -2 02clean_data/sample1_1_R2_paired_clean.fq.gz \

-p 12 -o 04salmon_quant/sample1_1

使用示例命令可以得到7个文件。quant.sf是根据transcript_id定量的结果，quant.genes.sf是根据gene_id定量的结果。

3、kallisto和salmon比较

对于免比对软件kallisto^[7]和salmon来说，比对和定量是一步完成的。对于salmon^[8]来说，可以通过添加-g参数，将原来的meta-feature(transcript_id)转换为我们想要的meta-feature(gene_id)，甚至他还可以自己判断测序类型。整体来说salmon表现更佳^[8]，而kallisto是science常用^[9]。

三、总结

上游分析包括质量检测、质量控制、比对、定量。质量检测和质量控制使用fastqc和Trimmomatic既可，比对和定量软件的选择比较多。比对软件STAR和HISAT2，定量软件featureCounts、RSEM和HTSeq-count。比对软件处理质控后的数据后可以得到sam和bam文件，其中记录到reads比对的位置、可变剪接、比对质量、配对reads比对到的位置等。定量软件基于比对软件给出的sam或bam文件进行定量。主要是根据比对到的位置进行定量，当然也会有一定的筛选条件。

最常用的组合套装是STAR-RSEM和HISAT2-featureCounts以及HISAT2-HTSeq-count。从整体上看，RSEM很全面的，因为它调用了STAR做联配，所以效率高速度快，而且这个组合输出的文件相当丰富，除了基于基因组和转录本比对的bam文件，其定量文件还包含count、TPM、FPKM等，为后续分析提供很大便利。相比之下HISAT2-featureCounts的定量结果只包含count和TPM，HISAT2-HTSeq-count的定量结果仅有count，比对结果只有sam，可以明显感觉到这两套组合的输出文件并不丰富，但胜在快捷、目的性强。如果我们不需要STAR-RSEM产生的诸多文件，那么这两套就足够了。

免比对软件kallisto和salmon对PC用户是首选，因为其可以跳过比对的步骤，直接拿质控后的数据进行定量，所以需要的内存和运算量比组合套装少得多。尤其是salmon可以自动判断测序类型，对于不是很了解整个流程和细节而又不想花时间学习的用户来说，使用salmon再好不过。

参考文献

[1] Sahraeian Sayed Mohammad Ebrahim,Mohiyuddin Marghoob,Sebra Robert,Tilgner Hagen,Afshar Pegah T,Au Kin Fai,Bani Asadi Narges,Gerstein Mark B,Wong Wing Hung,Snyder Michael P,Schadt Eric,Lam Hugo Y K. Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis.[J]. Nature communications,2017,8(1):59.

[2] 忘川水.RNA-seq转录组上游分析流程(2021年)[OL].https://zhuanlan.zhihu.com/p/369749492. 2021.

[3]Grüning Björn,Dale Ryan,Sjödin Andreas,Chapman Brad A,Rowe Jillian,Tomkins-Tinch Christopher H,Valieris Renan,Köster Johannes. Bioconda: sustainable and comprehensive software distribution for the life sciences.[J]. Nature methods,2018,15(7):475-476.

[4]刘永鑫Adam.数据的质量控制软件——fastQC[OL].https://blog.csdn.net/woodcorpse/article/details/106552332. 2018.

[5]Bolger, A. M., Lohse, M., & Usadel, B.Trimmomatic: A flexible trimmer for Illumina Sequence Data.[J]. Bioinformatics,2014,30(15):2114-20.

[6]Liao Y, Smyth GK and Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics,2014,30(7):923-30.

[7]Bray Nicolas L,Pimentel Harold,Melsted Páll,Pachter Lior. Near-optimal probabilistic RNA-seq quantification.[J]. Nature biotechnology,2016,34(5):525-7.

[8]Patro Rob,Duggal Geet,Love Michael I,Irizarry Rafael A,Kingsford Carl. Salmon provides fast and bias-aware quantification of transcript expression.[J]. Nature methods,2017,14(4):417-419.

[9]马省伟.使用salmon和kallisto进行RNA-seq定量[OL].https://blog.sciencenet.cn/blog-1094241-1133526.html. 2018.