12. 比对 定量

系统运维 2023-10-05 大猫 手机阅读
  1. 引入

      我们得到参考基因组以及测序序列后,接下来就是比对到参考基因组上.但是但常理说我们应该拿RNA的参考序列与测序数据比对,但是我们并没有RNA的参考序列,因此就只能拿测序数据比对到参考基因组上.实际还有一种做法就是用三代测序测完整的RNA序列,但是三代测序比较贵,因此不怎么考虑.

      使用的包是Hisat2Hisat2Hisat2.

  2. alignment

  接下来进入比对环节.我们进行参考基因组比对,需要建立IndexIndexIndex,上次我们使用的包是seqkitseqkitseqkit,这次我们使用hisat2hisat2hisat2.当然也可以使用STAR,只是这里使用了Hisat2

2.1 建立Index

  hisat2-build构建IndexIndexIndex,下面的代码是最简单的构建方式.这里如果基因组不在当前文件夹,可以使用ln指令,相当于设置一个快捷方式.语法:

ln [参数][源文件或目录][目标文件或目录]
    -s 软链接(符号链接)
ln -s ../data/genome.fa
#为../data/genome.fa在当前路径下设置一个软连接

hisat2-build [options]*  
#an example
hisat2-build genome.fa genome

  我们构建Index的时候,也可以提供给hisat2hisat2hisat2其他信息,以供IndexIndexIndex更好地构建.比如下面代码:

# 调用hisat2包的hisat2_extract_splice_sites.py从gtf(基因注释文件)中提取物种剪切位点的信息
hisat2_extract_splice_sites.py ../data/SWO_genes.gtf > genome.ss
1. 调用hisat2包的hisat2_extract_splice_sites.py从gtf(基因注释文件)中提取物种exon的信息
hisat2_extract_exons.py ../data/SWO_genes.gtf > genome.exon

  在这步后,我们重新构建Index.

hisat2-build --exon genome.exon --ss genome.ss ../data/genome.fa genome 1>histsat2-build.log 2>&1

  也可以写进脚本,然后让放进后台自动运行.

nohup sh step1.run_hist2.sh &

这里有一段关于hisat2_extract_exons的进一步解释Click.

2.2 比对

  我们可以先进行一个样本的比对.下面是代码,我们先解释参数意义.

hisat2 ‑‑new‑summary ‑p 2 ‑x genome ‑1 ../13.clean_data/CS1_1_1.fastq.gz ‑2 ../13.clean_data/CS1_1_2.fastq.gz ‑ S Cs1_1.sam ‑‑rna‑strandness RF

  • -x : 这里是参考基因组的索引名.对应于前面的genome.

  • -1 : 双端测序中左端readsreadsreads的fastq文件

  • -2 : 双端测序中右端readsreadsreads的fastq文件

  • -S : 要输出的比对结果文件路径,后缀为.sam

  • -p : 启动并行线程的数目

  • -new-summary : 以新样式打印对齐摘要，该样式对机器更友好

  • --rna‑strandness : 指定特定链的信息,默认为非链特异性文库.对于单末端reads,使用F或R."F"表示reads对应于转录本."R"是指reads对应于转录本的反向互补链.对于双末端reads,请使用FR或RF.

    • 使用此选项时,每个read比对结果都具有一个XS属性tag:"+"表示read属于基因组"+"链上的转录本."-"表示read属于基因组"-"链上的转录本.

  借此我想整理一下网上搜集的链特异性文库的tipstipstips.首先是正义链和反义链.

  • 正义链: 与模板链互补的那条链.这条链的碱基顺序和mRNA是一致的(除了其中的U/T替换).显然,这条非模板链对于mRNA是有意义的,存储着mRNA的编码信息,因此这条链也被称为编码链(coding strand)或者正义链(sense strand).
  • 反义链: mRNA转录时结合的那条链称之为模板链,也被称为非编码链和反义链.

反义链 = 非编码链 = 模板链 正义链 = 编码链 = 非模板链

  下文来自这位佬的定义,说得很清楚了.

  • 正链: 在参考基因组的fastqfastqfastq文件里,只给出了一条链的序列,这条链就是正链.这是硬性规定的,没有生物学意义.如果非要说有什么意义,那么可能就是为了世界各地的命名规则统一,使得研究更为方便.英文中通常将正负链写为forward/reverse strand,有时也会使用plus/minus strand表示.
  • 反链: 与正链互补的那条链.

------------------这里是我的个人理解--------------------

  我们在测序RNA时,会将RNA逆转录为DNA双链,这时我们将DNA分子拿去测序,得到的测序结果是不知道RNA由DNA双链的哪条链转录而来.这就无法区分cDNAcDNAcDNA中的正义链和反义链.此时就出现了链特异性文库.它的建库原理与普通转录组类似,不同之处在于合成第二链cDNA时,用dUTP代替dTTP,此时第二链cDNA上布满了含dUTP的位点.之后我们会用酶将第二链降解掉,此时就只剩下第一链.我们将第一链用于PCR扩增,然后进行测序,得到的reads(个人感觉这个reads是和第一链互补的)与参考基因组进行比对.如果reads比对上了,说明这条reads是位于正义链上.如果基因方向反方向就是反义链reads.

  诸如Ensembl和UCSC等数据库网站,它们更为关注的是基因的编码序列.因此当它们说一个基因位于正链上时,其含义是该基因的编码链位于正链上.也就是说该基因的RNA的模板链位于反链上.

  参考了这个人的博客Click以及Click.

----------------------有大佬纠正更好----------------------------

  回归正题,这个选项主要是标明文库是特异性文库还是非特异性文库.单端测序只需要标明R或F,如果是双端测序,就需要标明第一条链是转录本还是互补链.大多数公司都是使用RFRFRF.

  在指令运行完后,我们会将输出记录在日志文件里,打开日志文件内容如下:

  total pairs表示reads对数.9.09%的reads比对不上.2.9%的reads有多处不对(参考序列上可能有多处与reads相同的序列,这样就算不清楚该reads到底算哪个基因的表达量)。0.33%左右的reads比对上但是要遗弃(我们双端测序获得左右两端的reads,如果一端比对得上,另一端比对上的位置离这一端太远,比如5M以外,这就是不合理的情况,需要丢弃).我们用于定量的主要是比对上一次,且成对比对的(87.68%).   再看 Total unpaired reads,这是没有成对比上的数目.如果一端reads在染色体1上,另一端reads在染色体2上,这就是unpairedunpairedunpaired比对的情况.   总的比对率是94.45%.这算比较好的,如果是80%左右算比较差的.

image.png

  设置链特异性文库比非特异性文库比对率低是正常的.如果设置了,就只会按照设置的方向比对;如果没有设置,就会将readsreadsreads按正向比对一次,再反向比对一次,所以没有设置链特异性文库会更慢一点.

2.3 Sam文件

  这个运行时间会比较久,我们用后台运行,运行后可以得到23G的.sam.sam.sam文件.下文摘自wikiwikiwiki百科

SAM（Sequence Alignment Map，可直译为“序列比对地图”）是生物信息学中一种用于储存已比对到基因组上的序列信息的文件格式.

  SAM文件由两部分组成，头部区和主体区，都以tab分列.bam文件是sam文件的二进制格式,极其节约储存空间.

  转自Click.

2.3.1 头部区

  • @HD VN:1.0 SO:unsorted （排序类型） 头部区第一行:VN是格式版本;SO表示比对排序的类型,有unknown(default),unsorted，queryname和coordinate几种.对于coordinate，排序的主键是Alignments section的第三列"RNAME",其顺序由@SQ行的"SN"标识的顺序定义,次要排序键是Alignments section的第四列"POS"字段.对于RNAME和POS相等的比对,排列顺序则是任意的.

  • @SQ SN:contig1 LN:9401 (序列ID及长度) 参考序列名,这些参考序列决定了比对结果sort的顺序,SN是参考序列名;LN是参考序列长度;每个参考序列为一行. 例如：@SQ SN:NC_000067.6 LN:195471971

  • @RG ID:sample01 （样品基本信息） Read Group。1个sample的测序结果为1个Read Group；该sample可以有多个library的测序结果，可以利用bwa mem -R 加上去这些信息。 例如：@RG ID:ZX1_ID SM:ZX1 LB:PE400 PU:Illumina PL:Miseq ID：样品的ID号 SM：样品名 LB：文库名 PU：测序以 PL：测序平台 这些信息可以在形成sam文件时加入，ID是必须要有的后面是否添加看分析要求

  • @PG ID:bowtie2 PN:bowtie2 VN:2.0.0-beta7 （比对所使用的软件及版本） 例如： @PG ID:bwa PN:bwa VN:0.7.12-r1039 CL:bwa sampe -a 400 -f ZX1.sam -r @RG ID:ZX1_ID SM:ZX1 LB:PE400 PU:Illumina PL:Miseq …/0_Reference/Reference_Sequence.fa ZX_HQ_clean_R1.fq.sai ZX_HQ_clean_R2.fq.sai …/2_HQData/ZX_HQ_clean_R1.fq …/2_HQData/ZX_HQ_clean_R2.fq

这里的ID是bwa，PN是bwa，VN是0.7.12-r1039版本。CL可以认为是运行程序,@RG是上面RG表示的内容，后面是程序内容，这里的@GR内容是可以自己在运行程序是加入的

image.png

2.3.2 实体区

  该部分包含了11列必需字段,无效或者没有的字段一般用"0"或者"*"表示.建议直接看这个Click和这个.

主体部分有11个主列和1个可选列

  • QNAME 比对的序列名称 例如：M04650:84:000000000-B837R:1:1101:22699:1759（一条测序reads的名称）
  • FLAG Bwise FLAG（表明比对类型：paring，strand，mate strand等） 例如：99.我们可以用samtools flag flag值查看比对类型.
  • RENAME 比对上的参考序列名 例如：NC_000075.6
  • POS read比对到的参考序列"RNAME"最左侧的位置坐标,也是CIGAR中第一个比对标识"M"对应的最左侧碱基在参考序列的位置.未比对上的read在参考序列中没有坐标，此列标识为"0".
  • MAPQ 比对质量 例如：60
  • CIGAR read中每个碱基的比对情况 例如：87M
  • MRNM reads是成对的,这里是另外一条序列,比对上的参考序列名 例如：=
  • MPOS 1-Based leftmost Mate Position(相比于MRNM列来讲意思和POS差不多) 例如：124057667
  • ISIZE 插入片段长度 例如：200
  • SEQ 和参考序列在同一个链上比对的序列(若比对结果在负义链上,则序列是其反向重复序列，反向互补序列) 例如:ATTACTTGGCTGCT
  • QUAL 比对序列的质量（ASCII-33=Phred base quality）reads碱基质量值 例如:-8CCCGFCCCF7@E-
  • 可选的列 以TAG：TYPE：VALUE的形式提供额外的信息

image.png

2.4 将sam转为bam

  转换使用的工具是samtoolssamtoolssamtools.使用condacondaconda安装即可.这个软件包专门用于处理samsamsam文件.将samsamsam转为bambambam的指令称为samtools view.我们也可以用这个指令查看bam文件.

nohup samtools view -b -o CS1_1.bam -@ 4 Cs1_1.sam &

  我们说明一下那几个选项.

  • -b : 默认输出是 SAM 格式文件,该参数设置输出 BAM 格式
  • -o : 输出文件名
  • -@ : 进程个数
  1. 表达定量

  这里主要使用表达定量的三个软件htseq‑counthtseq‑counthtseq‑count和subreadsubreadsubread里的featureCountsfeatureCountsfeatureCounts.featureCountsfeatureCountsfeatureCounts有命令行版和RRR语言版.虽然看的教程使用的是RRR语言版,但是这里使用命令行版了.

3.1 featuresCounts的安装

  官网下看右边有直接的下载连接,我们下载Linux版的tar.gz,然后在自己路径下解压,插入环境变量即可.

image.png

3.2 featuresCounts的使用

  很明显要标注表达定量,需要比对结果(bambambam),以及基因组注释文件(gtfgtfgtf).这两个文件一起可以计算出每条基因上有多少个readsreadsreads.   看官方提示是最简单的方法,但是我好像点不进去官网文档.这里就一个个参数解释吧.

featureCounts -s 2 -T 4 -p --countReadPairs -M  -t exon -g gene_id  -a ../data/SWO_genes.gtf -o count.txt ../21.Mapping/CS1_1.bam

-s有关的看这篇博客Click

  • -s : 执行特定链的读取计数.应提供单个整数值(应用于所有输入文件)或一串逗号分隔值(应用于每个相应的输入文件).可能的值包括：0(非链式)、1(链式)和2(反向链式).默认值为0(即对所有输入文件执行非链式读取计数).这里是个人理解,绝大多数公司采用dUTP建库,测的RNA是反向链式
  • -T : 线程个数
  • -p : 如果指定,则文库会认为是含双端reads的文库.
  • --countReadPairs : 配合上面使用,只有当fragmentfragmentfragment(也就是两端一起)比对上了才会被计数,只有一个read会被忽略.
  • -M : 当一个fragmentfragmentfragment比对到多个地方也会计数.
  • -t : 指明feature−typefeature-typefeature−type,必须是gtf文件里有的featurefeaturefeature,当指明这个时,只有当reads在这些feature区域内才会被统计
  • -g : 在gtfgtfgtf文件中有id dentifierid dentifierid dentifier,当我们指明了id dentifierid dentifierid dentifier后,会将在id dentifierid dentifierid dentifier区域内的reads全部相加,作为id dentifierid dentifierid dentifier的统计结果.默认是gene_id,也就是说gene_id相同的reads会被相加到一起,这是基因层面的定量,我们也可以写transcript_id,这是对转录本层面的定量.

image.png

  这里我想阐明两个概念.feature指的是基因组区间的最小单位,比如exon;而metafeature可以看做是许多的feature构成的区间,比如属于同一个gene的外显子的组合.也可以说−g-g−g将featurefeaturefeature层面的结果统计到meta−featuremeta-featuremeta−feature层面.   之后,可以看到下面的图片,paired−endpaired-endpaired−end表明是双末端测序, count read pairs表示只计算双端一起比对上的.Multimapping reads表示当一个fragmentfragmentfragment比对到多个地方也会计数.

image.png

3.3 featureCounts 输出

  指令运行完成后,会生成两个文件一个是count.txtcount.txtcount.txt另一个是count.txt.summarycount.txt.summarycount.txt.summary.

3.3.1 summary

  我们先看count.txt.summarycount.txt.summarycount.txt.summary.   先来解释是什么意思. Assigned表示readsreadsreads比对到genesgenesgenes上的数目.unassigned表示没有分配到genes上的reads,下划线后接的是原因.unmappedunmappedunmapped表示根本没有比对上.NoFeaturesNoFeaturesNoFeatures表示比对到genes之间的序列.AmbiguityAmbiguityAmbiguity表示比对到多个地方,不知道归类给谁.

image.png

3.3.2 count.txt

  可以说输出文件很"工整".我们把列名认识一下.从左到右依次：

Geneid：基因的ensemble基因号；
Chr：多个外显子所在的染色体编号;
Start：多个外显子起始位点,与前面一一对应
End：多个外显子终止位点,与前面一一对应
Strand：正负链
Length：基因长度
sampleID：一列代表一个样本，数值表示比对到该基因上的read数目

image.png

  现在得到的表格内容只能称得上是raw matrix,我们还需要计算FPKM、TPM等,如下图.这个是老师写的,我自己没研究出来这么多小数位是怎么保存的....

image.png

  我自己写的代码如下,不过只计算了FPKMFPKMFPKM,需要再添加一列新列计算TPMTPMTPM.还在研究中.

sed "1,2d" count.txt | awk '{print $6 "t" $7 "t" "22706085"}' | awk '{print $1"t" $2"t" $3"t"  $2*1000000000/$1/$3}' | less -SN

image.png

3.4 merge

  上面的计算完成后,需要将多个样本的表达矩阵mergemergemerge起来.这里使用的是老师修改后的abundance_estimates_to_matrix.pl文件.走的代码如下:

perl ../software/abundance_estimates_to_matrix.pl --est_method featureCounts --quant_files genes.quant_files.txt --out_prefix gene

  之后生成如下文件:

image.png   我们可以打开一个看看,比如 genes.counts.matrix.可以看到行名是样本名、列名是基因编号.中间的值是reads比对上的次数counts.

image.png    TMM校正后生成的文件名带有TMM.下图是genes.TMM.TPM.matrix.这个文件也是经过两步标准化后,最终的表达量文件.

image.png

  1. 扩展

3.1 bam可视化

  要想可视化展示bambambam文件,可以使用软件IGV.我们直接去官网安装,可以选择无需安装Java和需要安装java的版本.   安装完成后,如下所示的界面.

image.png   要使用IGV,首先需要准备一些数据.基因组序列fa文件、基因组序列Index的fai文件、基因组注释文件gtf、比对结果bam文件、以及bam的Index文件.   为什么需要bam的Index文件呢?在原来的bambambam文件里,我们比对结果是按基因的名字排序的(或者没排序),但如果我们要可视化,就需要让结果在染色体上的比对位置排序(比如从1号染色体开始往后排2号染色体等).我们要使用的指令是 samtools sort.

nohup samtools sort -o CS1_1.sorted.bam CS1_1.bam 1>sort.log 2>&1 &

  排完序后,我们要对排序后的bam文件建立index.会获得CS1_1.sorted.bam.bai

samtools index CS1_1.sorted.bam

  我们用IGV打开genome.fagenome.fagenome.fa文件,然后选择染色体chr1.通过放大可以看到染色体的序列.

image.png

  接下来打开gtfgtfgtf文件.将范围缩小,可以看到蓝色条.这个蓝色条对应着一个参考基因组序列.

image.png

  点击蓝条,可以看到里面具体包含的类型.

image.png

  之后也可以打开bam文件,可以看到比对信息.中间灰色部分是比对的上的序列.

image.png   我们也可以使用IGV看覆盖深度.

image.png

image.png

3.2 不同软件的选择

  看下图,我们对转录组进行分析也可以使用其他软件,我们根据目的性来选择软件.比如stringtie适合 转录本重构,我们用测序的数据组装成转录本.hisat2 featureCount就是普通的转录组分析,序列比对、定量分析等.什么时候我们需要重构呢?比如我们测序的是长非编码RNA,它是不保守的,在不同品种、不同个体之间都有较大差异,因此用参考基因组比对效率不高.

  kallisto/salmon + DESeq2/edgeR没有比对就直接进行定量,这比hisat2 featureCount缺少一个质控步骤,比对率的高低是不清楚的,但是可以节省大量时间.

image.png

3.3 小节总结

image.png