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.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

  • 1
    点赞
  • 20
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值