这是RNA-Seq 上游分析的大致流程，比对+定量。当然实验目的若只需要定量已知基因，也可以选择free-alignment 的流程工具如kallisto/Salmon/Sailfish，其优点是可用于RNA-seq的基因表达的快速定量，但是对于小RNA和表达量低的基因分析效果并不好（2018年刚发表的一篇文章对free-alignment 的工具进行了质量评估，doi: >1 比对的是：使用idba_ud拼接的AER314-4raw_data基因组与转录组数据。
2 bowtie2做index（bowtie2使用conda安装）
建索引：bowtie2-build AER314-4_scaffoldfa AER314-4_scaffoldfa
3 reads mapping到参考基因组——tophat2软件：基于bowtie2(tophat安装见软件安装）
命令：tophat2 -p 12 -o AER314-4_output /home/test04/lyr/rna-seq/02align_out/AER314-4_scaffoldfa /home/test04/lyr/rna-seq/01data/YSH-qurRNA-42-314-4_L001_R1fastq /home/test04/lyr/rna-seq/01data/YSH-qurRNA-42-314-4_L001_R2fastq
4 然后就很顺利的跑出来结果了
使用公司服务器，12个线程，大概五个小时就跑完啦。
5 cufflink
[ Cufflinks输出结果
cufflinks的输入文件是sam或bam格式。并且sam或bam格式的文件必须排好序。（The SAM file supplied to Cufflinks must be sorted by reference position）Tophat的输出结果sam或bam已经排好了序。针对其他的未排序的sam或bam文件采用如下排序方式：
sort -k 3,3 -k 4,4n hitssam > hitssamsorted
1 transcriptsgtf
该文件包含Cufflinks的组装结果isoforms。前7列为标准的GTF格式，最后一列为attributes。其每一列的意义：
列数 列的名称 例子 描述
1 序列名 chrX 染色体或contig名; 2 来源 Cufflinks 产生该文件的程序名; 3 类型 exon 记录的类型，一般是transcript或exon; 4 起始 1 1-base的值; 5 结束 1000 结束位置; 6 得分 1000 ; 7 链 + Cufflinks猜测isoform来自参考序列的那一条链，一般是'+','-'或'';8 frame Cufflinks不去预测起始或终止密码子框的位置; 9 attributes 详见下
每一个GTF记录包含如下attributes：
Attribute 例子 描述
gene_idCUFF1Cufflinks的gene id;transcript_idCUFF11 Cufflinks的转录子 id; FPKM 101267 isoform水平上的丰度, F ragments P er K ilobase of exon model per M illion mapped fragments; frac 07647 保留着的一项，忽略即可，以后可能会取消这个;conf_lo 007 isoform丰度的95%置信区间的下边界，即 下边界值 = FPKM ( 10 - conf_lo );conf_hi 01102 isoform丰度的95%置信区间的上边界，即 上边界值 = FPKM ( 10 + conf_hi ); cov 100765 计算整个transcript上read的覆盖度;full_read_support yes 当使用 RABT assembly 时，该选项报告所有的introns和exons是否完全被reads所覆盖
2 ispformsfpkm_tracking
isoforms(可以理解为gene的各个外显子)的fpkm计算结果
3 genesfpkm_tracking
gene的fpkm计算结果Cuffmerge简介
Cuffmerge将各个Cufflinks生成的transcriptsgtf文件融合称为一个更加全面的transcripts注释结果文件mergedgtf。以利于用Cuffdiff来分析基因差异表达。
2 使用方法
$ cuffmerge [options]
输入文件为一个文本文件，是包含着GTF文件路径的list。常用例子：
$ cuffmerge -o /merged_asm -p 8 assembly_listtxt
3 使用参数
-h | --help
-o default: /merged_asm
将结果输出至该文件夹。
-g | --ref-gtf将该reference GTF一起融合到最终结果中。
-p | --num-threads defautl: 1
使用的CPU线程数
-s | --ref-sequence /该参数指向基因组DNA序列。如果是一个文件夹，则每个contig则是一个fasta文件；如果是一个fasta文件，则所有的contigs都需要在里面。Cuffmerge将使用该ref-sequence来帮助对transfrags分类，并排除repeats。比如transcripts包含一些小写碱基的将归类到repeats ]
cufflinks：
<1>命令：cufflinks -p 4 -o test_cuff /home/andengdi/lyr/rna-seq/02-align_out/test_output/accepted_hitsbam
流程及结果
5 用相同的方法将其他两个样本跑一遍。
subread是个套件,里面有subread aligner, subjunc aligner, featureCounts, exactSNP。subread aligner可以用于DNA-seq和RNA-seq；当用于RNA-seq时,subread只适用于差异分析；对于检测基因组变异如可变剪接之类的,需要reads的完全比对,这时候可以使用subjunc进行比对；在比对RNA-seq数据时,subread不会取检测exon-exon junctions的存在,只会把exon-spanning eads的最大可比对区域作为比对结果但是,如果只是进行差异分析的话,subread的结果足以进行subread的比对上reads可能会比subjunc多
HISAT2利用大量FM索引以覆盖整个基因组，使用改进的BWT算法，实现了更快的速度和更少的资源占用。
TopHat是一个基于Bowtie的RNA-Seq数据分析工具，它可以快速确认exon-exon剪切拼接事件。TopHat2比对有三个主要的过程，转录组的比对、基因组的比对、剪接比对，双端测序的Read首先每端数据要分别单独进行比对，然后考虑比对的片段长度以及方向将单独比对结果合并在一起形成双端比对结果。samtools是一个用于 *** 作sam和bam文件（通常是短序列比对工具如bwa,bowtie2,hisat2,tophat2等等产生的，具体格式可以在消息框输入“SAM”查看）的工具合集，包含有许多命令。以下是常用命令的介绍。
1View
view命令的主要功能是：将sam文件与bam文件互换；然后对bam文件进行各种 *** 作，比如数据的排序(sort)和提取(这些 *** 作 是对bam文件进行的，因而当输入为sam文件的时候，不能进行该 *** 作)；最后将排序或提取得到的数据输出为bam或sam（默认的）格式。
bam文件优点：bam文件为二进制文件，占用的磁盘空间比sam文本文件小；利用bam二进制文件的运算速度快。
view命令中，对sam文件头部（序列ID）的输入(-t或-T）和输出(-h)是单独的一些参数来控制的。
Usage: samtools view [options] <inbam>|<insam> [region1 []]
默认情况下不加 region，则是输出所有的 region
options:
-b output BAM
默认下输出是 SAM 格式文件，该参数设置输出 BAM 格式
-h print header for the SAM output
默认下输出的 sam 格式文件不带 header，该参数设定输出sam文件时带 header 信息
-H print header only (no alignments)
仅仅输出文件的头
-S input is SAM
默认下输入是 BAM 文件，若是输入是 SAM 文件，则最好加该参数，否则有时候会报错。
-u uncompressed BAM output (force -b)
该参数的使用需要有-b参数，能节约时间，但是需要更多磁盘空间。
-c Instead of printing the alignments, only count them and print the
total number All filter options, such as ‘-f’, ‘-F’ and ‘-q’ , are taken into account
过滤功统计功能
-c print only the count of matching records
-L FILE output alignments overlapping the input BED FILE [null]
-t FILE list of reference names and lengths (force -S) [null]
使用一个list文件来作为header的输入
-T FILE reference sequence file (force -S) [null]
使用序列fasta文件作为header的输入
-o FILE output file name [stdout]
-F INT filtering flag, 0 for unset [0]
Skip alignments with bits present in INT [0]
数字4代表该序列没有比对到参考序列上
数字8代表该序列的mate序列没有比对到参考序列上
过滤功能。如F12过滤只有双端map的
-q INT minimum mapping quality [0]
比对的最低质量值，一般认为20就为unique比对了，可以结合上述-bF参数使用使用提取特定的比对结果
例子：
将sam文件转换成bam文件
samtools view -bS abcsam > abcbam
BAM转换为SAM
samtools view -h -o outsam outbam
提取比对到参考序列上的比对结果
samtools view -bF 4 abcbam > abcFbam
提取paired reads中两条reads都比对到参考序列上的比对结果，只需要把两个4+8的值12作为过滤参数即可
samtools view -bF 12 abcbam > abcF12bam
提取没有比对到参考序列上的比对结果
samtools view -bf 4 abcbam > abcfbam
提取bam文件中比对到caffold1上的比对结果，并保存到sam文件格式
samtools view abcbam scaffold1 > scaffold1sam
提取scaffold1上能比对到30k到100k区域的比对结果
samtools view abcbam scaffold1:30000-100000 $gt; scaffold1_30k-100ksam
根据fasta文件，将 header 加入到 sam 或 bam 文件中
samtools view -T genomefasta -h scaffold1sam > scaffold1hsam
2Sort
sort对bam文件进行排序。一些软件需要sort的bam或者sam文件，如stringtie，所以必须要sort使用；求depth时，也必须要sort；
Usage: samtools sort [-n] [-m <maxMem>] <inbam> <outprefix>
-m 内存参数默认下是 500,000,000 即500M（不支持K，M，G等缩写）。对于处理大数据时，如果内存够用，则设置大点的值，以节约时间。
-n 设定排序方式按short reads的ID排序。默认下是按序列在fasta文件中的顺序（即header）和序列从左往右的位点排序。
例子：
samtools sort acceptbam acceptsort最终产生acceptsortbam
3merge
将2个或2个以上的已经sort了的bam文件融合成一个bam文件。融合后的文件已经sort过了的。
Usage: samtools merge [-nr] [-h inhsam] <outbam> <in1bam> <in2bam>[]
Options: -n sort by read names
-r attach RG tag (inferred from file names)
-u uncompressed BAM output
-f overwrite the output BAM if exist
-1 compress level 1
-R STR merge file in the specified region STR [all]
-h FILE copy the header in FILE to <outbam> [in1bam]
例子：
4index
必须对bam文件进行默认情况下的排序后，才能进行index。否则会报错。
建立索引后将产生后缀为bai的文件，用于快速的随机处理。很多情况下需要有bai文件的存在，特别是显示序列比对情况下。比如samtool的tview命令就需要；gbrowse2显示reads的比对图形的时候也需要。IGV显示比对情况也需要。
Usage: samtools index <inbam> [outindex]
例子：
以下两种命令结果一样
$ samtools index abcsortbam
$ samtools index abcsortbam abcsortbambai
5faidx
对fasta文件建立索引,生成的索引文件以fai后缀结尾。该命令也能依据索引文件快速提取fasta文件中的某一条（子）序列
Usage: samtools faidx <inbam> [ []]
对基因组文件建立索引，方便提取序列。
例子：$ samtools faidx genomefasta
由于有索引文件，可以使用以下命令很快从基因组中提取到fasta格式的子序列
$ samtools faidx genomefasta scffold_10 > scaffold_10fasta
6tview
tview能直观的显示出reads比对基因组的情况，和基因组浏览器有点类似。
需要事先利用利用上面讲的sort和建index命令执行完后，用下述命令。
Usage: samtools tview <alnbam> [reffasta]
出参考基因组的时候，会在第一排显示参考基因组的序列，否则，第一排全用N表示。
按下 g ，则提示输入要到达基因组的某一个位点。例子“scaffold_10:1000"表示到达第
10号scaffold的第1000个碱基位点处。
使用H(左）J（上）K（下）L（右）移动显示界面。大写字母移动快，小写字母移动慢。
使用空格建向左快速移动（和 L 类似），使用Backspace键向左快速移动（和 H 类似）。
Ctrl+H 向左移动1kb碱基距离； Ctrl+L 向右移动1kb碱基距离
可以用颜色标注比对质量，碱基质量，核苷酸等。30～40的碱基质量或比对质量使用白色表示；
20～30；10～20绿色；0～10蓝色。
使用点号''切换显示碱基和点号；使用r切换显示read name等
还有很多其它的使用说明，具体按 ？ 键来查看。
7flagstat
给出BAM文件的比对结果
Usage: samtools flagstat <inbam>
$ samtools flagstat examplebam
11945742 + 0 in total (QC-passed reads + QC-failed reads)
#总共的reads数
0 + 0 duplicates
7536364 + 0 mapped (6309%:-nan%)
#总体上reads的匹配率
11945742 + 0 paired in sequencing
#有多少reads是属于paired reads
5972871 + 0 read1
#reads1中的reads数
5972871 + 0 read2
#reads2中的reads数
6412042 + 0 properly paired (5368%:-nan%)
#完美匹配的reads数：比对到同一条参考序列，并且两条reads之间的距离符合设置的阈值
6899708 + 0 with itself and mate mapped
#paired reads中两条都比对到参考序列上的reads数
636656 + 0 singletons (533%:-nan%)
#单独一条匹配到参考序列上的reads数，和上一个相加，则是总的匹配上的reads数。
469868 + 0 with mate mapped to a different chr
#paired reads中两条分别比对到两条不同的参考序列的reads数
243047 + 0 with mate mapped to a different chr (mapQ>=5)
#同上一个，只是其中比对质量>=5的reads的数量
8depth
得到每个碱基位点的测序深度,并输出到标准输出，所以要用大于号追加到一个文件。
Usage: bam2depth [-r reg] [-q baseQthres] [-Q mapQthres] [-b inbed] <in1bam> []
-r 后面跟染色体号（region）
-q ：计算深度时要求测序碱基质量最低质量值
-Q ：计算深度时要求比对的最低质量值
注意：做depth之前必须做samtools index；
例子
samtools depth acceptbam >depth
9其他命令
reheader：替换bam文件的头
$ samtools reheader <inheadersam> <inbam>
idxstats ：统计一个表格，4列，分别为”序列名，序列长度，比对上的reads数，unmapped reads number。第4列应该是paired reads中有一端能匹配到该scaffold上，而另外一端不匹配到任何scaffolds上的reads数。
$ samtools idxstats <alnbam>
rmdup：将由PCR duplicates获得的reads去掉，并只保留最高比对质量的read。
Usage: samtools rmdup [-sS]
-s 对single-end reads。默认情况下，只对paired-end reads
-S 将Paired-end reads作为single-end reads处理。
10 将bam文件转换为fastq文件
有时候，我们需要提取出比对到一段参考序列的reads，进行小范围的分析，以利于debug等。这时需要将bam或sam文件转换为fastq格式。
该网站提供了一个bam转换为fastq的程序：>
宅在家两个多月，不知不觉已经是春天了，也许距离返校的日子更近了吧
变异 ，指的是实际测序数据与国际规定的参考基因组之间的区别。很多变异其实只是造成人类多样性的原因。 突变 ，指的是那些与疾病相关的变异。
举个例子：ENSEMBL等规定的人类参考基因组文件某位置是AAAAA，然后一个人实际测序得到的序列为AGCAA，那么相比于参考基因组，这个人就有2个变异位点。对于第2个位置，如果查看所有已知的测序，绝大部分人都是G，说明是参考基因组出现了问题，这个变异就不能称作突变。对于第3个位置，如果查看所有已知的测序，绝大部分人都是A，而恰好有一个人不是A，但他是个患者，那么这个变异就是突变了。
SNP（single nucleotide polymorphism）：单核苷酸多态性。 个体间基因组DNA序列同一位置单个核苷酸变异（替换、插入或缺失）所引起的多态性。在人类基因组中SNP分布普遍并且密度较大，总数超过107， 平均每300bp（也有说1kbp）就有一个SNP。或称单核苷酸位点变异SNV。
INDEL（insertion-deletion）：插入和缺失。 基因组上小片段（>50bp）的插入或缺失。
CNV（copy number variation）：基因组拷贝数变异。 基因组中大片段的DNA形成非正常的拷贝数量。比如一个基因在染色体的一条染色单体上的数目为1，但是在染色体复制过程中，复制结束后该基因在染色单体数目由1变成了2或者n。它发生的频率远远高于染色体结构变异，并且整个基因组中覆盖的核苷酸总数大大超过SNP的总数。
SV（structure variation）：结构变异。 染色体大片段的插入与缺失，染色体内部的某区域发生翻转颠换，两条染色体之间发生重组。
一般情况下只分析SNP，其它类型的变异分析有难度或不准确。
来自两个不同个体的DNA片段AAGCCTA和AAGCTTA为等位基因。几乎所有常见的SNP位点只有两个等位基因。
在人体中，SNP的发生机率大约是01%，也就是每1000个碱基对就可能有一个SNP（密度高）。对疾病发生和药物治疗有重大影响的SNP，估计只占数以百万计SNP的很小一部分。
SNP位点的分布是不均匀的，在非转录序列比在转录序列更常见。编码区的单核苷酸多态性——编码 SNP（coding SNP，cSNP）也有同义和非同义两种类型，非同义SNP会改变蛋白质的氨基酸序列。基因非编码区、基因间隔区的SNP仍然可能影响转录因子结合、剪接等过程。
从演化的观点来看，SNP具有相当程度的稳定性，即使经过代代相传，SNP所引起的改变却不大，因此可用以研究族群演化。
HISAT2 是一款利用改进的BWT算法进行序列比对的软件。由约翰霍普金斯大学计算生物学中心（CCB at JHU）开发，是TopHat的升级版本，速度提高了50倍。利用 HISAT2 + StringTie 流程，可以快速地分析转录组测序数据，获得每个基因和转录本的表达量。
首先需要构建参考基因组索引用于下一步的比对。HISAT2提供了两个脚本用于从基因组注释GTF文件中提取剪接位点和外显子位置，基于这些特征，可以使 RNA-Seq reads 比对更加准确。然后再进行reads mapping。
比对结果：
SAM（sequence Alignment/mapping）数据格式是目前高通量测序中存放比对数据的标准格式。BAM是SAM的二进制格式。使用samtools将sam文件转化为bam文件，并进行排序。
SAM文件：
vcf格式（Variant Call Format）是存储变异位点的标准格式，用于记录variants（SNP / InDel）。BCF是VCF的二进制文件。
stats统计文件：
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）