RNA-Seq数据分析，质控，STAR比对步骤

一、需要下载的软件

fastqc、multiqc、star
GitHub上STAR的网站：https://github.com/alexdobin/STAR
STAR 我下的软件是2.6.1a版本的，自己根据需要下载合适的版本。

  wget -c https://github.com/alexdobin/STAR/archive/refs/tags/2.6.1a.tar.gztar -xzvf 2.6.1a.tar.gz#解压chmod -R 755 STAR-2.6.1a#给执行权限

二、数据下载

公司返回的数据需要下载到自己服务器上的路径。

三、md5check

##md5check
cd /mnt/raid62/BetaCoV/Person/zengwanqin/data/DBW_RNA_seq/X101SC21030878-Z03-F004/2.cleandata
ls > name.txt
vim name.txt #删除非样本名
for i in `cat /mnt/raid62/BetaCoV/Person/zengwanqin/data/DBW_RNA_seq/X101SC21030878-Z03-F004/2.cleandata/name.txt`
do
cd $i
md5sum -c MD5* >> /mnt/raid62/BetaCoV/Person/zengwanqin/data/DBW_RNA_seq/X101SC21030878-Z03-F004/2.cleandata/md5/md5_check_result.txt
cd ..
done

这一步后会得到一个md5_check_result.txt文本文件，内容如下：

全是OK就说明无误，可进行后续分析。

四、fastqc、multiqc

cd /mnt/raid62/BetaCoV/Person/zengwanqin/data/DBW_RNA_seq/X101SC21030878-Z03-F004/2.cleandata
for i in `cat name.txt`
do
bash /mnt/raid61/Personal_data/lizhidan/scripts/rnaPipeline/bash/01_01_fastqc.bash \
-i /mnt/raid62/BetaCoV/Person/zengwanqin/data/DBW_RNA_seq/X101SC21030878-Z03-F004/2.cleandata/$i \
-o /mnt/data2/zengwanqin/data/fastq/dbw_fastqc/ -p 8 -n $i
donecd /mnt/data2/zengwanqin/data/fastq/dbw_fastqc
multiqc .

可得到一个multiqc_report.html网页版，打开可看质控结果。

1、General Statistics：所有样本数据基本情况统计

%Dups——重复reads的比例
%GC——GC含量占总碱基的比例，比例越小越好
Length——测序长度
M Seqs——总测序量（单位：millions）

2、 Sequence Quality Histograms：每个read各位置碱基的平均测序质量

横坐标——碱基的位置
纵坐标——质量分数
质量分数=-10log10p（p代表错误率），所以当质量分数为40的时候，p就是0.0001。此时说明测序质量非常好。
绿色区间——质量很好，
橙色区间——质量合理。
红色区间——质量不好。

3、PerSequence Quality Scores 具有平均质量分数的reads的数量

横坐标——平均序列质量分数
纵坐标——reads数
绿色区间——质量很好
橙色区间——质量合理
红色区间——质量不好
当峰值小于27时——warning
当峰值小于20时——fail

4、Per Base Sequence Content ：每个read各位置碱基ATCG的比列

对所有reads的每一个位置，统计ATCG四种碱基的分布。

横坐标——碱基位置
纵坐标——样本
%T——红色
%C——蓝色
%A——绿色
%G——紫色

reads每个位置的颜色显示由4种颜色的比例混合而成，哪一个碱基的比例大，则趋近于这个碱基所代表的颜色。

正常情况下每个位置每种碱基出现的概率是相近的。
如果ATGC在任何位置的差值大于10%——warning
如果ATGC在任何位置的差值大于20%——fail

5、Per Sequence GC Content ：reads的平均GC含量

横坐标——GC含量百分比
纵坐标——数量

正常的样本的GC含量曲线会趋近于正态分布曲线，曲线形状的偏差往往是由于文库的污染或是部分reads构成的子集有偏差（overrepresented reads）。形状接近正态但偏离理论分布的情况提示我们可能有系统偏差。

偏离理论分布的reads超过15%时——warning
偏离理论分布的reads超过30%时——fail

6 、Per Base N Content ：每条reads各位置N碱基含量比例

当测序仪器不能辨别某条reads的某个位置到底是什么碱基时，就会产生“N”，统计N的比率。正常情况下，N值非常小。

横坐标——read中的位置
纵坐标——N的数量比
当任意位置的N的比例超过5%——warning
当任意位置的N的比例超过20%——fail

7 、Sequence Duplication Levels：每个序列的相对重复水平

横坐标：每个序列的相对重复水平
纵坐标：在文库中的比例
当非unique的reads占总数的比例大于20%时——warning
当非unique的reads占总数的比例大于50%时——fail

测序深度越高，越容易产生一定程度的duplication，这是正常的现象，但如果duplication的程度很高，就提示我们可能有bias的存在。

8、Overrepresented sequences：文库中过表达序列的比例

横坐标——过表达序列的比例
纵坐标——样本
过表达序列的比例>0.1%——warning
过表达序列的比例>1%——warning

一条序列的重复数，因为一个转录组中有非常多的转录本，一条序列再怎么多也不太会占整个转录组的一小部分（比如1%），如果出现这种情况，不是这种转录本巨量表达，就是样品被污染。这个模块列出来大于全部转录组1%的reads序列，但是因为用的是前100,000条reads，所以其实参考意义不大。

9 、Adapter Content 接头含量

横坐标——碱基位置
纵坐标——占序列的百分比
- 大于5%——warning
- 大于10%——fail

10、Status Checks 总结版

五、STAR

cd /mnt/raid62/BetaCoV/Person/zengwanqin/data/DBW_RNA_seq/X101SC21030878-Z03-F004/2.cleandata
cat name.txt|xargs -P 4 -I {} STAR --runThreadN 4 \--genomeDir /mnt/raid64/ref_genomes/MusMus/release93/STAR_index_2.7.1a \#构建好的index--outSAMtype BAM SortedByCoordinate \--outBAMcompression 9 \--readFilesIn /mnt/raid62/BetaCoV/Person/zengwanqin/data/DBW_RNA_seq/X101SC21030878-Z03-F004/2.cleandata/{}/{}_1.clean.fq.gz  /mnt/raid62/BetaCoV/Person/zengwanqin/data/DBW_RNA_seq/X101SC21030878-Z03-F004/2.cleandata/{}/{}_2.clean.fq.gz \--readFilesCommand zcat \--outFileNamePrefix /mnt/data2/zengwanqin/data/dbw_align/{}. \--outSJfilterOverhangMin 12 12 12 12 \--alignSJoverhangMin 3 \--alignSJDBoverhangMin 3 \--chimSegmentMin 12 \--chimScoreMin 2 \--chimScoreSeparation 10 \--chimJunctionOverhangMin 12 \--outFilterMultimapNmax 1 \ --chimOutType Junctions SeparateSAMold

参数说明

--runThreadN 设置线程数
--runMode alignReads : 默认就是比对模式，可以不填写
--genomeDir: 索引文件夹
--readFilesIn FASTA/Q文件路径
--readFilesCommand zcat: 如果输入格式是gz结尾，那么需要加上zcat, 否则会报错
--outSAMtype: 输出SAM文件的格式，是否排序
--outBAMsortingThreadN: SAM排序成BAM时调用线程数
几个额外要说明的点：- 由于物种的组装的复杂性，存在一些为组装上的片段，这些片段不需要放在参考序列中，尤其是可变单倍型(alternative haplotypes)
- 如果基因组的contig过多，超过5000，你需要用 --genomeChrBinNbits=min(18,log2[max(GenomeLength/NumberOfReferences,ReadLength)]) 降低RAM消耗
- 选择最新的注释文件，人类和小鼠常在http://www.gencodegenes.org下载，植物的可信基因组见http://plants.ensembl.org
- 如果没有设置--sjdbGTFfile或--sjdbFileChrStartEnd，就不需要设置--sjdbOverhang

step1、建立基因组索引文件index
– runThreadN NumberOfThreads
–runMode genomeGenerate
–genomeDir /path/to/genomeDir #index输出的路径
–genomeFastaFiles /path/to/genome/fasta1 #参考基因组序列，解压缩文件
–sjdbGTFfile /path/to/annotations.gtf #参考基因组注释文件
–sjdbOverhang ReadLength-1 #这个是reads长度的最大值减1，默认是100
（其他参数像产生gff3或者可变剪接比对等，需要加特殊参数）step2、mapping
–runThreadN NumberOfThreads #线程
–genomeDir /path/to/genomeDir #index目录
–readFilesIn /path/to/read1 [/path/to/read2 ] #paired reads文件
–outFileNamePrefix /path/to/output/dir/prefix #默认结果输出pwd（–outReadsUnmapped 因为需要去除一些假性重复的序列，例如rRNA元件。所以需要比对到 RepBase_human_database_file ）

这一步可得到以下文件，可用于后续分析。

"SJ.out.tab"存放的高可信的剪切位点，每一列的含义如下

第一列: 染色体
第二列: 内含子起始（以1为基）
第三列: 内含子结束（以1为基）
第四列：所在链，1(+)，2(-)
第五列: 内含子类型，0表示不是下面的任何一种功能，1表示GT/AG, 2表示:GT/AC,3表示GC/AG,4表示GT/GC,5表示AT/GC,6表示GT/AT
第六列: 是否是已知的注释
第七列: 有多少唯一联配支持
第八列: 有多少多重联配支持
第九列: maximum spliced alignment overhang, 这个比较难以翻译，指的是当短读比对到剪切位点时，中间会被分开，另一边能和基因组匹配的数目，例如ACGTACGT----------ACGT，就是4或者8，取决于方向。
控制过滤的参数为–outSJfilter*系列，其中–outSJfilterCountUniqueMin 3 1 1 1表示4类内含子唯一匹配的read支持数至少为3,1,1,1, 而–outSJfilterCountTotalMin 3 1 1 1则表示4类内含子唯一匹配和多重匹配read的支持数和，至少为3,1,1,1。如果你设置的–outSJfilterReads Unique，那么上面两者是等价的，当然默认情况下是All

参考网页：
https://www.jianshu.com/p/85da4dcc6020
https://www.jianshu.com/p/fa388b21d1de