一、需要下载的软件
- 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.gz
tar -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
done
cd /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