转录组学习四（参考基因组及gtf注释探究）

转录组学习一（软件安装）转录组学习二（数据下载）转录组学习三（数据质控）转录组学习四（参考基因组及gtf注释探究）转录组学习五（reads的比对与samtools排序）转录组学习六（reads计数与标准化）转录组学习七（差异基因分析）转录组学习八（功能富集分析）任务在UCSC下载hg19参考基因组，从gencode数据库下载基因注释文件，并且用IGV去查看你感兴趣的基因的结构，比如TP53,KRAS,EGFR等等。下载ENSEMBL，NCBI的gtf，也导入IGV看看，截图基因结构。了解IGV常识。<font color=orange>参考基因组的选择</font><font size=3>首先</font>转录组分析流程里十分重要的一步就是对reads进行mapping（比对） 到参考基因组或者参考转录组进而定量进行后续分析。对于非模式物种来说就是用二代三代测序数据组装成参考转录组。而对于有参考基因组的物种来说，就十分方便直接从数据库中下载参考基因组即可。由于此次作业数据是为对人类和小鼠进行的rna-seq数据。故直接下载人类和小鼠的参考全基因组数据即可。人类的参考基因组：所了解的主要有<code><font color=red>NCBI版、ENSEMBL版、UCSC版</font></code>。只有各自格式上的不同，序列则是相同的。根据在生信菜鸟团我的基因组（五）:测试数据及参考基因组的准备对参考基因组的选择上描述：这个对新手来说，是一个很大的坑，hg19、GRCH37、 ensembl 75这3种基因组版本应该是大家见得比较多的了，国际通用的人类参考基因组，其实他们储存的是同样的fasta序列，只是分别对应着三种国际生物信息学数据库资源收集存储单位，即NCBI，UCSC及ENSEMBL各自发布的基因组信息而已。有一些参考基因组比较小众，存储的序列也不一样，比如BGI做的炎黄基因组，还有DNA双螺旋结构提出者沃森（Watson）的基因组，还有2016年发表在nature上面的号称最完善的韩国人做的基因组。前期我们先不考虑这些小众基因组，主要就下载hg19和hg38，都是UCSC提供的，虽然hg38相比hg19来说，做了很多改进，优点也不少，但因为目前为止很多注释信息都是针对于hg19的坐标系统来的，我们就都下载了，正好自己探究一下。也顺便下载一个小鼠的最新版参考基因组吧，反正比对也就是睡个觉的功夫，顺便分析一下结果，看看比对率是不是很低。选择下载UCSC的hg38，UCSC网址(https://genome.ucsc.edu/index.html)mkdir Database && cd Databasewget http://10.10.0.8/hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gzgunzip hg38.fa.gz参考基因组的组装形式：分别是primary, toplevel。三种重复序列处理方式：分别是unmasked(DNA)、soft-masked(dna_sm)和masked(dna_rm)。一般<font color = red >选择dna.primary或dna_sm.primary。</font>根据参考基因组和基因注释文件对这些文件的选择：为什么选择PrimaryPrimary assembly contains all toplevel sequence regions excluding haplotypes and patches. This file is best used for performing sequence similarity searches where patch and haplotype sequences would confuse analysis.为什么不选择maskedMasked基因组是指所有重复区和低复杂区被N代替的基因组序列，比对时就不会有reads比对到这些区域。一般不推荐用masked的基因组，因为它造成了信息的丢失，由此带来的一个问题是uniquely比对到masked基因组上的reads实际上可能不是unique的。而且masked基因组还会带来比对错误，使得在允许错配的情况下，本来来自重复区的reads比对到基因组的其它位置。 另外检测重复区和低复杂区的软件不可能是完美的，这就造成遮盖住的重复序列和低复杂区并不一定是100%准确和敏感的。soft-masked基因组是指把所有重复区和低复杂区的序列用小写字母标出的基因组，由于主要的比对软件，比如BWA、bowtie2等都忽略这些soft-mask，直接把小写字母当做大写字母比对，所以使用soft-masked基因组的比对效果和使用unmasked基因组的比对效果是相同的。<font color=orange>基因注释文件的解释</font>两个重要的文件：GTF文件和GFF文件对于GTF文件来说：其基本情况如图：

image1 havana gene 11869 14409 . + . gene_id "ENSG00000223972"; gene_version "5"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; havana_gene "OTTHUMG00000000961"; havana_gene_version "2";1 havana transcript 11869 14409 . + . gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; havana_gene "OTTHUMG00000000961"; havana_gene_version "2"; transcript_name "DDX11L1-002"; transcript_source "havana"; transcript_biotype "processed_transcript"; havana_transcript "OTTHUMT00000362751"; havana_transcript_version "1"; tag "basic"; transcript_support_level "1";1 havana exon 11869 12227 . + . gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "1"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; havana_gene "OTTHUMG00000000961"; havana_gene_version "2"; transcript_name "DDX11L1-002"; transcript_source "havana"; transcript_biotype "processed_transcript"; havana_transcript "OTTHUMT00000362751"; havana_transcript_version "1"; exon_id "ENSE00002234944"; exon_version "1"; tag "basic"; transcript_support_level "1";1 havana exon 12613 12721 . + . gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; havana_gene "OTTHUMG00000000961"; havana_gene_version "2"; 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"; transcript_support_level "1";1 havana exon 13221 14409 . + . gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "3"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; havana_gene "OTTHUMG00000000961"; havana_gene_version "2"; transcript_name "DDX11L1-002"; transcript_source "havana"; transcript_biotype "processed_transcript"; havana_transcript "OTTHUMT00000362751"; havana_transcript_version "1"; exon_id "ENSE00002312635"; exon_version "1"; tag "basic"; transcript_support_level "1";1 havana transcript 12010 13670 . + . gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; havana_gene "OTTHUMG00000000961"; havana_gene_version "2"; transcript_name "DDX11L1-001"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; havana_transcript "OTTHUMT00000002844"; havana_transcript_version "2"; tag "basic"; transcript_support_level "NA";1 havana exon 12010 12057 . + . gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "1"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; havana_gene "OTTHUMG00000000961"; havana_gene_version "2"; transcript_name "DDX11L1-001"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; havana_transcript "OTTHUMT00000002844"; havana_transcript_version "2"; exon_id "ENSE00001948541"; exon_version "1"; tag "basic"; transcript_support_level "NA";1. seq_id：序列的编号，一般为chr或者scanfold编号2. source: 注释的来源，一般为数据库或者注释的机构，如果未知，则用点“.”代替3. type:注释信息的类型，比如Gene、cDNA、mRNA、CDS等;4. start: 该基因或转录本在参考序列上的起始位置；(从1开始，包含);5. end: 该基因或转录本在参考序列上的终止位置；(从1开始，包含);6. score: 得分，数字，是注释信息可能性的说明，可以是序列相似性比对时的E-values值或者基因预测是的P-values值，.表示为空;7. strand: 该基因或转录本位于参考序列的正链(+)或负链(-)上;8. phase: 仅对注释类型为“CDS”有效，表示起始编码的位置，有效值为0、12. (对于编码蛋白质的CDS来说，本列指定下一个密码子开始的位置。每3个核苷酸翻译一个氨基酸，从0开始，CDS的起始位置，除以3，余数就是这个值，，表示到达下一个密码子需要跳过的碱基个数。该编码区第一个密码子的位置，取值0,1,2。0表示该编码框的第一个密码子第一个碱基位于其5’末端；1表示该编码框的第一个密码子的第一个碱基位于该编码区外；2表示该编码框的第一个密码子的第一、二个碱基位于该编码区外；如果Feature为CDS时，必须指明具体值。)；9. attributes: 一个包含众多属性的列表，格式为“标签＝值”(tag=value)，以多个键值对组成的注释信息描述，键与值之间用“=”，不同的键值用“；”隔开，一个键可以有多个值，不同值用“,”分割。注意如果描述中包括tab键以及“，= ；”，要用URL转义规则进行转义，如tab键用 代替。键是区分大小写的，以大写字母开头的键是预先定义好的，在后面可能被其他注释信息所调用。ID：注释信息的编号，在一个GFF文件中必须唯一；name：注释信息的名称，可以重复；Alias：别名；Parent > >Indicates：该注释所属的注释，值为注释信息的编号，比如外显子所属的转录组编号，转录组所属的基因的编号。Parent指明feature所从属的上一级ID，用于将exons聚集成transcript，将transripts聚集成gene，值可以为多个；Target 指明比对的目标区域，一般用于表明序列的比对结果。格式为 “target_idstart end [strand] ,其中strand是可选的(“+”或”-”),target_id中如果包含空格，则要转换成’ '。Gap：T比对结果的gap信息，和Target一起，用于表明序列的比对结果。Derives_from：Note:备注;Dbxref:数据库索引。<font color=orange>基因注释文件的探究</font>统计每条染色体上基因的数目、(CDS、transcript、等也类似)。==PERL==open IN,"$ARGV[0]" or die "$!";while (<IN>){ chomp; next if(/^#/); @line=split; if($line[2] eq "gene"){ @chr=$line[0]; ###注意，此@列表是可变的。分别为【1】、【2】..【100】 foreach (@chr){ $count{$_}+=1; } }}foreach $chr (sort keys %count){ print "$chr\t=>\t$count{$chr}\n";}grep与AWKzcat Homo_sapiens.GRCh38.87.chr.gtf.gz | grep -v "^#" | awk '$3~/gene/{print $0}' |less -S##此是显示所有包含gene的行。其中$0是所有列。 grep -v "^#" Homo_sapiens.GRCh38.87.chr.gtf| awk '$3~/gene/{print $0}' | awk '$0~/protein_coding/ {print $1,"\t",$3}' |sort | uniq -c |less -S ###统计 gene_biotype为protein_coding的数目，运用到了awk{sum+=$2}END{print sum}此awkgrep -v "^#" Homo_sapiens.GRCh38.87.chr.gtf | awk '$3~/gene/{print $0}' | awk '$0~/gene_biotype "protein_coding"/{print$1,"\t",$3}'| uniq -c | awk '{sum+=$1}END{print sum}'统计物种GTF文件里gene_biotype的类型及各自的数量==PERL==#!/usr/bin/perl -wwhile(<>){ next if (/^#/); @line=split; if($line[2] eq "gene"){ if (/gene_biotype "(.*?)";/){ push @biotype,$1; } }}foreach (@biotype){ $count{$_}+=1;}foreach(sort keys %count){ print"$_\t=>\t$count{$_}\n";}==perl单行==zcat Homo_sapiens.GRCh38.87.chr.gtf.gz |perl -alne '{next unless $F[2] eq "gene" ;/gene_biotype "(.*?)";/;print $1}' |sort |uniq -c###?最小匹配数结果：

image