附录 E — 基础问题 21—25
本组保留原题和原稿参考答案。题库目前收录 1—25 题；历史仪器参数和软件用法仍需逐项更新。
E.1 BBQ100-21
问题描述
Hello大家好！我们今天又见面了！
今天是第21题，我们接着之前的题目，继续学习与SAM/BAM有关的内容。今天要学习的内容是SAM/BAM文件的附加信息。
1. 基础导引部分
我们先给大家举个例子，这是一个human的全基因组测序比对的SAM文件的11列以后的信息。第11列之前学习过了是reads的质量值，那么后面的若干标记比如MD:Z:145等等这些符号是什么意思呢？
图1 SAM文件的11列以后的信息截图
我把上面图中的部分行的信息放到这里，供大家查阅（11列以后的内容要一直向右拖拽)
ST-E00126:128:HJFLHCCXX:2:1206:8105:9730 99 chr1 11670 1 145M = 11898 315 AGGTGAAGCCCTGGAGATTCTTATTAGTGATTTGGGCTGGGGCCTGGCCATGTGTATTTTTTTAAATTTCCACTGATGATTTTGCTGCATGGCCGGTGTTGAGAATGACTGCGCAAATTTGCCGGATTTCCTTTGCTGTTCCTGC KKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKFKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKFKKFKKKKKKKKKKKKKKKKKKKFFKKKKKKKKKKFKKKKKKKKKKKKKKKFAK MD:Z:145 PG:Z:MarkDuplicates XG:i:0 NM:i:0 XM:i:0 XN:i:0 XO:i:0 AS:i:0 XS:i:0 YS:i:0 YT:Z:CP
ST-E00126:128:HJFLHCCXX:2:2107:22820:18520 99 chr1 11682 1 145M = 11920 325 GGAGATTCTTATTAGTGATTTCGGCTGGTGCCTGGCCATGTGTATTTTTTTAAATTTCCACTGATGATTTTGCTGCATGGCCGGTGTTGAGAATGACTGCGCAAATTTGCCGGATTTCCTTTGCTGTTCCTGCATGTAGTTTAAA KKKKAAKKAFFKKKKKKFKFKKKFKKKKKKKKKKFKFKKKKKKKKKKKKKKFKFFKKKKKKFAAKAKKKKKKKKKKKKFFKKKFFFKKFKFFKKKKKKKKFFFFFKKKKKKK7<FFKKKKKKAFK<F<<7<AA,,7AA<7F7AA< MD:Z:21G6G116 PG:Z:MarkDuplicates XG:i:0 NM:i:2 XM:i:2 XN:i:0 XO:i:0 AS:i:-12 XS:i:-12 YS:i:-6 YT:Z:CP
ST-E00126:128:HJFLHCCXX:2:1210:9110:60026 163 chr1 11703 1 87M = 11840 282 GGGCTGGGGCCTGGCCATGTGTATTTTTTTAAATTTCCACTGATGATTTTGCTGCATGGCCGGTGTTGAGAATGACTGTGCAAATTT 7FKAFKFKFFFKKF<FKKKFKKKFKK7F<KFFKFKKKKKKKFFFF,FKKKKKKKFFKKKKK(7<,AAK<F7AAFKKFKFKFF<A<7< MD:Z:78C8 PG:Z:MarkDuplicates XG:i:0 NM:i:1 XM:i:1 XN:i:0 XO:i:0 AS:i:-5 XS:i:-5 YS:i:-33 YT:Z:CP
ST-E00126:128:HJFLHCCXX:2:2101:7425:68324 99 chr1 11708 1 145M = 11923 302 GGGGCCTGGCCATGTGTATTTTTTTAAATTTCCACTGATGATTTTGCTGCATGGCCGGTGTTGAGAATGACTGCGCAAATTTGCCGGATTTCCTTTGCTGTTCCTGCATGTAGTTTAAACGAGATTGCCAGCACCGGGTATCATT KKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK<FKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKAKKKKK<FKKKKFKKKKKKKKK7<AKFFAFKFF<KKKKKFKK<FK<7F,AFKFFA MD:Z:145 PG:Z:MarkDuplicates XG:i:0 NM:i:0 XM:i:0 XN:i:0 XO:i:0 AS:i:0 XS:i:0 YS:i:0 YT:Z:CP
ST-E00126:128:HJFLHCCXX:2:2210:15382:54752 163 chr1 11714 1 87M = 11866 297 TGGCCATGTGTATTTTTTTAAATTTCCACTGATGATTTTGCTGCATGGCCGGTGTTGAGAATGACTGCGCAAATTTGCCGGATTTCC KKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKFKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK MD:Z:87 PG:Z:MarkDuplicates XG:i:0 NM:i:0 XM:i:0 XN:i:0 XO:i:0 AS:i:0 XS:i:0 YS:i:-5 YT:Z:CP 一般呢，我们都把11列以后的内容称为可选择区域（optional fields），这个区域所有的格式都必须是TAG:TYPE:VALUE的形式，比如MD:Z:145就是一个符合规范的可选区域的值。
根据SAM格式官方文档的信息，我们需要记住以下内容：
1. 所有的TAG都是2个字母，一般情况下都是大写字母。并且TAG在1行的比对结果中只能出现1次。
2. 所有的TYPE都是单字母，大小写敏感，它是用来定义后面VALUE的类型；
3. VALUE可长可短，但是需要和之前的TYPE相呼应。
关于TYPE不同字母对应的不同数据类型，把SAM的官方文档贴一下，共大家参考。其中，最常用的就是i（带符号的数字）；Z（可直接输出字符串，可以包含空格）；
图2 TYPE的字母与不同数据类型之间的对应关系
2.常用的TAG
那么常用的TAG都有哪些，都代表什么含义呢？要知道，不同的比对软件可能会在SAM文件的后面加上不同的TAG，所以我们在查询TAG含义的时候一定要从所用比对软件的官方文档中去查找。而SAM文件的header部分又包含了@PG字符段可以帮助我们还原比对软件的参数设置，因此我们拿到一个SAM文件就可以通过查阅文档的方式了解TAG的基本信息。
使用samtools可以查看sam文件的header部分
samtools view -H test.sam比如，我们这里的@PG内容如下
@PG ID:bowtie2-5DEB9F7A PN:bowtie2 VN:2.2.5 CL:"/home/biotools/bowtie2-2.2.5/bowtie2-align-s --wrapper basic-0 -p 4 --phred33 -x /lustre/user/reference/hg19/hg19_combine -S ./tmp.data/fastq/genome-sequence.sam -1 ./tmp.data/fastq/genome-sequence_L3_1_trim5.fastq -2 ./tmp.data/fastq/genome-sequence_L3_2_trim5_92.fastq"3. 提问环节
我们今天的问题很简单，请根据bowtie2的官方文档，解释下面的比对信息：
ST-E00126:128:HJFLHCCXX:2:2107:22820:18520 99 chr1 11682 1 145M = 11920 325 GGAGATTCTTATTAGTGATTTCGGCTGGTGCCTGGCCATGTGTATTTTTTTAAATTTCCACTGATGATTTTGCTGCATGGCCGGTGTTGAGAATGACTGCGCAAATTTGCCGGATTTCCTTTGCTGTTCCTGCATGTAGTTTAAA KKKKAAKKAFFKKKKKKFKFKKKFKKKKKKKKKKFKFKKKKKKKKKKKKKKFKFFKKKKKKFAAKAKKKKKKKKKKKKFFKKKFFFKKFKFFKKKKKKKKFFFFFKKKKKKK7<FFKKKKKKAFK<F<<7<AA,,7AA<7F7AA< MD:Z:21G6G116 XG:i:0 NM:i:2 XM:i:2 XN:i:0 XO:i:0 AS:i:-12 XS:i:-12 YS:i:-6 YT:Z:CP
Anwser
1. ST-E00126:128:HJFLHCCXX:2:2107:22820:18520
> 序列名称，比对片段的编号，通常包括测序平台的信息
2. 99
> Flag值
3. chr1
> 回帖到的染色体名称
4. 11682
> 比对到染色体上的具体位置（比对到正链最左边bp的位置点）
5. 1
> 比对的质量值，叫做MAPQ，MAPQ=-10 * log10{mapping出错的概率}
6. 145M
> CIGAR值，描述具体的比对情况
7. =
> pair reads中与该序列配对的read所mapping到的参考序列，如果没有mapping到同一条参考序列上，则用“*”代替。
8. 11920
> pair reads中与该序列配对的read所mapping到的参考序列的具体位置
9. 325
> 通过分析pair reads mapping到同一条参考序列上位置的推断得到fragment的长度
10. GGAGA....TTAAA
> read序列信息
11. KKKKA....F7AA<
> read序列测序每一bp的质量值
12. MD:Z:21G6G116
> MD:Z:表示在比对过程中有mismatch的情况，后面字符串表示mismatch的具体位置
13. XG:i:0
> XG:i有gap的存在，后面数字表示gap的总长度（read和reference上的都计算在内）
14. NM:i:2
> 编辑距离，为了将read map到reference上，对read进行单核苷酸编辑（替换、插入以及删除）的最小长度
15. XM:i:2
> mismatche的具体数目
16. XN:i:0
> 序列覆盖区的参考基因组上不确定的base数
17. XO:i:0
> gap的具体数目
18. AS:i:-12
> 比对分数，允许负值，局部比对最终可以大于0，但是全局比对中不会
19. XS:i:-12
> 比对过程中出现的比最终报告分数（AS:i:-12）高的比对值，同样允许负值，局部比对最终可以大于0，但是全局比对中不会。当一条序列能够同时比对到多个位点，且出现连续局部相似度极高的情况下会出现这种情况。
20. YS:i:-6
> 与该序列配对的pair read的比对分数
21. YT:Z:CP
> YT:Z:代表pair-read的比对情况，“UU”代表没有配对的read； "CP"代表序列为pair reads之一，pair align cordantly；"DP"表序列为pair reads之一，pair align discordantly；"UP"代表序列为pair reads之一，但是pair没有比对到参考基因组上。
参考资料：
E.2 BBQ100-22
问题描述
Hello 大家好！
前面的若干问题，我们一直在围绕着SAM文件的记录格式做了详细地讨论，我相信大家通过我们的问题，跟随我们学习的思路已经掌握了SAM文件作为标准的比对格式的合理性以及相关特点。
1. 背景介绍和数据下载
SAM文件不但记录了reads详细的mapping信息，还记录了reads的原始信息，内容很是全面。这样很好，但也存在很多问题：
比如我的原始FASTQ文件是100G，那么我的SAM文件一定是大于100G的，也就是占用了更多空间； mapping的结果是没有排序，无论是按reads的name排序还是按在基因组上的位置排序，都没有。所以默认的SAM输出文件是乱序的，处理很不方便； mapping的结果不能进行随机访问，什么是随机访问呢？举个例子就是说对于一个SAM文件我不能快速地访问比如chr1 10000 - 200000这个区域的所有reads的mapping情况。 基于以上这3个问题，BAM文件就出现了，并且完美解决了上面3个问题。为了方便我们今天的展示和说明，我为大家准备了1个很小的SAM文件，大约只有4MB，请大家下载下来并完成我们的相关问题。
SAM测试文件的baidu盘下载地址
链接：https://pan.baidu.com/s/15gVVYPRu3VbF_uKbJUUGrA 密码：2drn 同时，我们今天要使用的工具是Linux下的samtools，请没有Linux的老铁去安装Linux（我们马上就会有教程出来）；请没有安装samtools软件的老铁使用conda安装需要软件，教程可以移步（用Anaconda快速搭建生物信息学分析平台）
2. 思路讲解
BAM文件是SAM文件的一种压缩格式，也是最常用的一种比对结果的压缩格式。它一般可以将SAM文件压缩到只有原来的20~30%大小，并且使用非常方便。
同时，对于BAM文件，我们一般还会进行排序，根据不同的需要，我们排序的方法一般有2种：第1种是按照mapping到的参考基因组的坐标上下游顺序来排序，是samtools的默认排序方法；第2种是按照reads name进行排序，需要增加一个-n参数。
对于一个已经排序好的BAM文件，我们通常会建立索引文件，后缀名一般是在BAM文件名的后面多个“.bai“。有了BAM以及索引文件的出现，我们就可以随机访问任意一段染色体区域的BAM文件。
3. 提出问题
那么我们今天的问题也很简单，就是使用samtools工具对我们的测试数据test.sam文件进行操作。具体要求如下：
1. 使用samtools view 命令查看test.sam的header，请记录各条染色体的长度；同时告知这个test.sam文件是使用哪种mapping软件进行mapping的？
查看header中的@PG ID，显示使用的mapping软件是bowtie2，header中显示各条染色体的长度如下图：
答1. 各条染色体的长度
2. 使用samtools view命令将test.sam文件转换成test.bam文件，并保留header区域，写出命令并记录test.sam，test.bam的文件大小。
使用的命令如下：( -b：输出文件为bam格式; -h,输出中包含header信息)
samtools view -b -h test.sam > test.bam
结果显示test.sam文件大小是3.9M，而test.bam文件则是660K，bam文件会小很多； 3. 使用less命令分别查看test.sam，test.bam文件，为什么bam文件会输出乱码？使用samtools view命令再试试看？
less命令可以正常查看sam文件，但是不能正常查看bam文件，因为bam文件是二进制文件，所以需要使用
samtools view test.bam来查看bam文件。4. 使用samtools sort命令对test.bam文件进行排序，输出文件名为test_sort.bam，并记录文件大小。 首先解释一下samtools sort命令：
sort命令的使用：
samtools sort [-l level] [-m maxMem] [-o out.bam] [-O format] [-n] [-T tmpprefix] [-@ threads]
[in.sam|in.bam]
参数：
-l INT 设置输出文件压缩等级。0-9，0是不压缩，9是压缩等级最高。不设置此参数时，使用默认
压缩等级；
-m INT 设置每个线程运行时的内存大小，可以使用K，M和G表示内存大小。
-n 设置按照read名称进行排序；
-o FILE 设置最终排序后的输出文件名；
-T PREFIX 设置临时文件的前缀；
-O FORMAT 设置最终输出的文件格式，可以是bam，sam或者cram，默认为bam；
-@ INT 设置排序和压缩是的线程数量，默认是单线程。对test.bam进行排序，不压缩，默认线程，设置最终输出名称为test_sort.bam，
默认临时文件前缀，默认输出bam文件，默认单线程；
使用命令如下：
samtools sort -o test_sort.bam test.bam
结果显示 test.bam 660K; test_sort.bam 660K 也就是排序之后的bam文件大小不变。5. 使用samtools index 对test_sort.bam建立index，写出命令并记录其文件大小。
samtools index [-bc] [-m INT] <in.bam> [out.index]
参数：
-b 创建bai索引文献(默认);
-c 创建csi索引文献;
-m INT 创建csi索引文献，最小间隔值2^INT;
> samtools index test_sort.bam
> ls -hs
结果： test_sort.bam.bai 4.0K6. 使用samtools tview使用下面的命令查看chr1:160000-160100区域的比对情况，并截图
使用的命令如下：
samtools tview -p chr1:160000-160100 test_sort.bam 
答2. chr1:160000-160100区域的比对情况
4. 参考资料
资料1：本次主要是对samtools的一个应用，我建议大家直接看samtools的说明文档，比如对于view功能，直接在命令行敲击samtools view，再按回车就能出现说明文档，如下图所示
图1 samtools view的说明文档
5. 多说几句话
大家以后要用的软件种类非常多，不可能所有的软件你都学过，总有一个从不会到会的过程。在学习使用各种软件的过程中要学会类比，要学会推理，要想清楚我们的input是什么output是什么。不能乱搞一气，想不明白其中的道理，就像一个黑盒子，最后就是你不知道扔进去是什么，也不知道扔出来是什么，这不完蛋了？
至于软件使用方法的学习，一定要多看官方的说明文档！入门的时候看看别人的介绍或者是指导资料什么的尚可，但是一定有了一定基础以后一定要多阅读官方的说明文档，受益无穷的！
E.3 BBQ100-23
问题描述
Hello大家好！我们今天又见面了！
我们通过前期的22个问题，从数据的简单质控，到测序数据的mapping，再到mapping后的SAM文件都有了一个比较清楚的认识。那么说了半天的mapping问题，一直都是在以DNA进行举例，RNA的比对我们都还没有谈。那么今天我们就来简单谈谈RNA序列的mapping，尤其是真核生物的RNA序列比对。
1. RNA与DNA结构的不同
一般来说，DNA的mapping比较容易，因为DNA在基因上是连续的，直接回贴到基因组就可以找到相应的定位。就比如我们常用的Whole Genome Sequence（WGS）即全基因组测序；或者是我们所说的ChIP-Seq即染色体免疫共沉淀测序都是直接对DNA进行建库测序，其测序结果都是FASTQ文件，直接用bowtie2，bwa比对到基因组就可以拿到标准的SAM文件。
但是RNA就不一样了，真核生物的RNA需要经过复杂的加工过程。在细胞中RNA层面的调控至少可以分成2个大的阶段co-transcription（转录的同时） 和 post-transcription（转录以后）其中的调控机制也有很多。
对我们mapping影响最大的因素是：真核生物转录出来的初步的mRNA都是带有intron（内含子）的，随后都需要在co-transcription（转录的同时） 或post-transcription（转录以后）阶段通过：1. alternative splicing（可变剪切）剪切掉intron；2.polyA尾巴； 3.加5’的帽子结构。这3个步骤，将不成熟的mRNA变为最终成熟的mRNA再转运出核，行使功能。
图1 通过可变剪切同1个基因可以形成多种蛋白（https://en.wikipedia.org/wiki/Alternative_splicing）
2.RNA比对的常用软件
目前大家最常用的转录组比对软件有下面几个：
- tophat2，应用最广泛的比对软件，但是速度很慢，已经基本被淘汰了，大约需要4~5G内存就能运行；
- hisat2，tophat2的原班人马搞得新一代转录组比对软件，比对速度大大提高，我强烈推荐，大约需要4~5G内存就能运行；
- STAR，非常适合于大量数据的并行计算，速度非常快，对于同时有参考基因组和参考转录组的物种，比对的准确率很高，不过index很大，至少需要30G以上内存才能运行。
3.提出问题
问题1：如果你有一套标准的polyA捕获得到的RNA-Seq测序数据，对reads进行了前处理工作与质量控制工作，但是你的比对策略为：先尝试mapping，把能mapping到基因组上的reads都先mapping；然后把不能进行mapping的reads进行一定规则的拆分，再进行第二轮mapping，从而解决跨intron区域的问题（以上为tophat的mapping策略）。请问，这样mapping的最大问题是什么？（提示，需要知道一些假基因的概念！）
首先解释一下假基因，假基因（Pseudogenes）是指是一类染色体上的基因片段。假基因的序列通常与
对应的基因相似，但至少是丧失了 一部分功能，基因不能表达或其编码的蛋白质没有功能；这种基因
在基因组上的分布非常普遍，那么在假基因普遍存在的情况下上述比对策略就会受到假基因的干扰，会有
很多基因比对到假基因上。问题2：在human中，是不是所有的蛋白基因（protein coding gene）都含有intron？
并不是，SRY基因是人体Y染色体上的一段基因，该基因是决定男性睾丸发育的主要基因，存在于Y染色体
的短臂末端上，该基因只有一个exon。
答2 SRY基因结构,没有intron.
问题3：在human中，是不是所有的蛋白基因的成熟mRNA都有polyA尾巴？
并不是，组蛋白mRNA末端就没有polyA尾巴。E.4 BBQ100-24
问题描述
Hello大家好！我们又见面了！
在第23问的时候，我们开始学习了转录组相关的生物信息学。在转录组分析的过程中，往往需要基因注释文件，通常是GFF或者是GTF文件，那么这两个文件的内容是什么？有什么特点呢？这是我们今天要探索的问题。
1. 我们为什么需要基因注释文件？
图1. 通过对外显子（exon）的可变剪切，同1个基因可以形成多种蛋白（https://en.wikipedia.org/wiki/Alternative_splicing）
我们的gene在基因组上的结构不是连续的，而是exon-intron-exon（exon=外显子，intron=内含子）分隔开的。基因要表达，首先会先发生转录过程，转录出包含intron的pre-mRNA序列，然后再进行可变剪切，加5’帽子，3’ PolyA尾巴等一系列复杂的加工过程才会形成成熟的mRNA。
在进行转录组序列比对，尤其是mRNA序列比对的时候，经常需要处理跨越两个exon之间的reads，所以在进行序列比对的时候往往需要对基因组有一个注释，告诉比对软件哪个位置是gene的exon，哪个地方是gene的intron，这个就是我们所说的基因注释信息。所以，它里面核心内容就是一大堆gene在基因组上的坐标，以及这个gene本身的一些属性。
2. GTF与GFF都是基因注释文件
GTF = General Transfer Format
GFF = General Feature Format
GFF有若干个版本，简单来说，GTF是GFF文件的其中一个版本，我们一般认为GTF文件就是GFF 2.0版本的内容。一个标准的GTF/GFF2.0文件需要包括9列内容，一个简单的示意图如下：
图2. 1个标准的GTF格式文件，文件不包括前面的行号
# 所有的列必须用TAB分隔，总共有9列内容，第9列是补充列；
# 补充列的内容可以为空，但是前面8列必须有内容，如果想表达空的概念，则需要用"."；
# 第1列 seqname
染色体的名称，需要与genome FASTA文件中的染色体名对应，别一个用"chr1"一个用"Chr1"；
# 第2列 source
注释来自哪里，比如图2表示来自NCBI RefSeq数据库；
# 第3列 feature
此行的注释类型，一般有exon，CDS，stop_codon, start_codon等等；
# 第4，5列 start，end
此行注释的起始和终止位置，标准的GTF/GFF都是以1为染色体的起点（1-based system）;
注意！无论这个gene是正链还是负链，start的坐标都小于end坐标；
# 第6列 score
一般存放打分值，比如拼装的可信度之类。下载的官方注释文件一般为0.0
# 第7列 strand
正链基因标记为 "+", 负链基因标记为 "-";
# 第8列 frame
只可能是0,1,2这3个值,表示与CDS中codon的相对位置；
0表示，这个region的第1bp就是正好是codon 三连密码子的第1个碱基；
1表示，这个region的第2bp就是正好是codon 三连密码子的第1个碱基；
2表示，这个region的第3bp就是正好是codon 三连密码子的第1个碱基；
# 第9列 attribute
一般会记录 gene_id 与transcript_id;
这一列是可选列，可以增加很多内容。在程序处理过程中，相同的attribute会合并在一起处理。
比如，所有gene_id=SGIP1的行都会先汇总在一起，表示1个基因。3.提出问题
既然是这样，我们今天就思考2个小问题，1个比较偏理论，1个是比较具体。
- 你认为GTF/GFF的文件格式设计合理吗？为什么？
并不是非常合理，这种格式虽然包括了注释需要的全部信息，但是同一个基因不同的elements并不在
一行，在mapping的时候需要用循环一行一行去判断该element是否还同属一个基因，比较耗费和内存，
如果不选择GTF文件，而是选择下图所示的“all files from selected table”文件，这种文件的注释信息的组合方式与GTF不同，是以一个基因为一行，包括了这个基因中的各个elements，这样的注释方法使比对过程更加便捷。 
答1. 另一种注释文件
all files from selected table 文件内容示例：
#bin name chrom strand txStart txEnd cdsStart cdsEnd exonCount exonStarts exonEnds score name2 cdsStartStat cdsEndStat exonFrames
1251 NM_004261.4 chr1 - 87328127 87380048 87329156 87379794 5 87328127,87333735,87346344,87368963,87379710, 87329288,87333785,87346408,87369131,87380048, 0 SELENOF cmpl cmpl 0,1,0,0,0,- 如果告知，transcript_id 为NM001308203.1，gene_id 为SGIP1, 在转录本上的坐标为101，那么对应基因组的坐标是多少？请写出答案与简要程序思路。注释信息如下：
chr1 hg19_ncbiRefSeq exon 66999252 66999355 0.000000 + . gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq start_codon 67000042 67000044 0.000000 + . gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq CDS 67000042 67000051 0.000000 + 0 gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq exon 66999929 67000051 0.000000 + . gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq CDS 67091530 67091593 0.000000 + 2 gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq exon 67091530 67091593 0.000000 + . gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq CDS 67098753 67098777 0.000000 + 1 gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq exon 67098753 67098777 0.000000 + . gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq CDS 67105460 67105516 0.000000 + 0 gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq exon 67105460 67105516 0.000000 + . gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq CDS 67108493 67108547 0.000000 + 0 gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq exon 67108493 67108547 0.000000 + . gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq CDS 67109227 67109402 0.000000 + 2 gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq exon 67109227 67109402 0.000000 + . gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq CDS 67136678 67136702 0.000000 + 0 gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq exon 67136678 67136702 0.000000 + . gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq CDS 67137627 67137678 0.000000 + 2 gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq exon 67137627 67137678 0.000000 + . gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq CDS 67138964 67139049 0.000000 + 1 gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq exon 67138964 67139049 0.000000 + . gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq CDS 67142687 67142779 0.000000 + 2 gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq exon 67142687 67142779 0.000000 + . gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq CDS 67145361 67145435 0.000000 + 2 gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq exon 67145361 67145435 0.000000 + . gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq CDS 67154831 67154958 0.000000 + 2 gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq exon 67154831 67154958 0.000000 + . gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq CDS 67155873 67155999 0.000000 + 0 gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq exon 67155873 67155999 0.000000 + . gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq CDS 67160122 67160187 0.000000 + 2 gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq exon 67160122 67160187 0.000000 + . gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq CDS 67184977 67185088 0.000000 + 2 gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq exon 67184977 67185088 0.000000 + . gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq CDS 67194947 67195102 0.000000 + 1 gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq exon 67194947 67195102 0.000000 + . gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq CDS 67199431 67199563 0.000000 + 1 gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq exon 67199431 67199563 0.000000 + . gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq CDS 67205018 67205220 0.000000 + 0 gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq exon 67205018 67205220 0.000000 + . gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq CDS 67206341 67206405 0.000000 + 1 gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq exon 67206341 67206405 0.000000 + . gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq CDS 67206955 67207119 0.000000 + 2 gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq exon 67206955 67207119 0.000000 + . gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq CDS 67208756 67208775 0.000000 + 2 gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq stop_codon 67208776 67208778 0.000000 + . gene_id "SGIP1"; transcript_id "NM_001308203.1";
chr1 hg19_ncbiRefSeq exon 67208756 67216822 0.000000判断的思路是根据注释信息计算chr1上的exon的分布与长度，转录本上的坐标值代表了该片段之前的exon长度，那么以该题为例，chr1上该基因第一个Exon的长度是 6699355-66999252+1=104， 转录本坐标是101，说明该片段就是map到这个exon上的，坐标是6699355+101-1=6699455
E.5 BBQ100-25
问题描述
Hello everyone! 我们又见面了！
我们昨天刚刚回答了GTF/GFF文件到底是什么的问题，那么我们今天尝试回答GTF/GFF文件是怎么来的，从哪里能够下载。
根据BBQ24问的介绍，我们知道了GTF是GFF的2.0版本。（以后我们把这两者到简称为GTF文件好了！）主要记录的内容是基因组上基因的结构，包括哪里是exon，哪里是intron，哪里是CDS，哪里是start codon，哪里是stop codon. 在提问的环节，我们也设置了问题让大家去吐槽GTF文件糟糕的的文件结构（希望大家能吐得开心！）。今天我们来为大家介绍GTF文件是怎么来的。
1. 什么是参考GTF/GFF文件
针对一些已经通过测序计划，拼装好基因组序列信息的物种，一般会同时提供其转录组的注释信息，也就是我们所说的GTF/GFF文件。常用的模式生物，比如human（人），mouse（小鼠），rat（大鼠），chicken（鸡），lizard（蜥蜴），Arabidopsis thaliana（拟南芥）等等都已经有非常好的 全基因组参考序列（FASTA文件），以及转录组注释信息（GTF或GFF文件）。因此，直接到能够提供下载地址的网站上下载就好了。图1给大家展示了已经公布参考基因组哺乳动物的系统发生树，大家可以看看人类和哪种动物演化距离最近。
图1 已经公布参考基因组的哺乳动物系统发生树（http://asia.ensembl.org/info/about/speciestree.html）
多说一句，这个下载下来的FASTA文件就是我们所谓的参考基因组，需要用这个文件去构建mapping的index；下载下来的GTF文件是转录组注释信息，一般在计算表达量的时候需要提供。
2.什么是拼装转录本
在进行转录组分析的过程中，我们经常会听到一句话叫“用XX软件拼装转录本”，这句话是什么意思呢？不都有参考转录组了，还要拼个啥转录组？
第1个方面，对于无参考基因组，无参考转录组的物种来说，往往需要通过RNA-Seq的数据自己拼出参考转录组，然后再进行下游的数据分析。所以，对于无参分析来说，往往需要自己拼装转录本，生成自己的参考转录组信息，也包括注释信息（GTF文件）。
第2个方面，对于有非常好注释的基因组，例如human来说。不同的细胞条件可能不同，一些永生化的细胞系往往都具有“癌症”的特征，这些细胞的转录组，基因组或多或少都有结构的变异，以及转录本的差异。
举个例子，有一个geneA，在参考转录组中注释的是从chr1:100015000进行转录，但是在另外的细胞系中有可能就是从chr1:98015050进行转录。这种转录起始和终点的不同是很常见的现象。因此，有时候为了比较严谨地进行下游序列分析，是需要根据已知的参考转录组以及测序数据对其进行一个修饰。常用的软件是cufflinks，stringtie等等。
不过，对于有参考转录组的物种，一般情况下，我还是建议不要去自己拼转录本，也不要去做什么所谓的修正，意义不大。除非你研究的体系非常特殊。
3. 从哪里下载参考转录组GTF/GFF文件呢？
关于下载GTF/GFF文件的内容，给大家介绍2个最常用的网站：
1个是UCSC genome browser （UCSC Genome Browser-网址链接） 1个是Ensembl（Ensembl 网址链接） 注意：对于植物的Ensembl网站（Ensembl Plants）
3.1 从UCSC genome browser下载human的GTF文件
1. 打开UCSC genome browser网站 （图3.1-1）
2. 在Tools里选择 Table Browser(图3.1-2）
3. 打开Table Browser以后，设置相关的需要内容（图3.1-3）
4. 点击get output即可下载
# hg19 = human genome 19是常用的human参考基因组版本号；
# RefSeq gene是全部经过人工检查过的gene注释文件；
图3.1-1 打开UCSC genome browser网站
| 25 图3 |
|---|
| 25 图4 |
| 图3.1-3 打开Table Browser以后，设置相关的需要内容 |
| *3.2 Ensembl下载human的GTF文件** 下载之前我必须跟大家提个醒。 |
对于动物相关的信息都请访问Ensembl的动物站：http://www.ensembl.org/index.html 对于植物相关的信息都请访问Ensembl的植物站：http://plants.ensembl.org/index.html |
| 们在这里还是以下载human hg19版本的GTF文件为例，操作步骤如下： |
. 登陆Ensembl网站，并跳转到hg19版本界面 （图3.2-1） . 继续选择跳转到hg19版本界面（图3.2-2） . 在hg19版本的Ensembl界面中选择download（图3.2-3） . 在download页面中选择Download a sequence or region （图3.2-4） . 在左边栏选择 FTP download 然后选择下载 GTF文件（图3.2-5） . 选择注释好的GTF进行下载（图3.2-6） |
| 25 图5 |
| 图3.2-1 登陆Ensembl网站，并跳转到hg19版本界面 |
图3.2-2 继续选择跳转到hg19版本界面
 |
| 图3.2-3 在hg19版本的Ensembl界面中选择download |
图3.2-5 在左边栏选择 FTP download 然后选择下载 GTF文件
图3.2-6 选择注释好的GTF进行下载
4. 提问环节
1. 请按照文中教程分别从UCSC Genome Browser，以及Ensembl网站上下载hg19的转录组注释的GTF格式文件。
下载的文件及压缩文件大小如下：
2. 下载这两个文件解压缩以后的大小是否有差异，差异大不大？
两个网站下载的GTF文件大小差异较大，解压之后hg19_RefSeq_GTF_UCSC文件126M，Homo_sapiens.GRCh37.87.chr.gtf文件大小是1.2G。3. 解压并使用Linux less命令打开这两个文件，观察这两个文件的transcript_id以及gene_id是否相同，再找找看有哪些其他地方的不同。
两个文件的transcript_id以及gene_id均不相同，Ensembl网站下载的gtf文件transcript_id以及
gene_id均以ENSG和ENSG开头，全称是Ensembl Transcript ID和Ensembl Gene ID，而UCSC网站
下载的gtf文件transcript_id以及gene_id均以NM开头,在NCBI数据库中代表mRNA；除此之外，UCSC
的注释信息非常简练，而Ensemble网站中的gtf文件注释信息相比于UCSC更加全面，但也存在冗余，NCBI
数据库中的基因注释被验证的比例更大，所以在比对策略上可以选择先使用UCSC的GTF文件筛选目标基因，
再利用Ensemble数据库的GTF文件找到更加详细的注释。Ensembl gtf文件：
1 havana exon 12613 12721 . + . gene_id "ENSG00000223972"; gene_version "4"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "2"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene"; transcript_name "DDX11L1-002"; transcript_source "havana"; transcript_biotype "processed_transcript"; havana_transcript "OTTHUMT00000362751"; havana_transcript_version "1"; exon_id "ENSE00003582793"; exon_version "1"; tag "basic";UCSC gtf文件：
chr1 hg19_ncbiRefSeq exon 66999929 67000051 0.000000 + . gene_id "NM_001308203.1"; transcript_id "NM_001308203.1";