gffread 不仅可以实现GTF与GFF的互相转换,而且还可以对GFF文件进行过滤处理。可以直接读取GTF文件。
相关database:最好在ucsc上下载。或者cellranger的database包,
#下载gtf/gff文档及hg19文件
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_19/gencode.v19.annotation.gff3.gz
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz
#使用gunzip及tar分别解压,并合并下载内容获得hg19.fa
cat *.fa >hg19.fa
用法:
Usage:
gffread <input_gff> [-g <genomic_seqs_fasta> | <dir>][-s <seq_info.fsize>]
[-o <outfile.gff>] [-t <tname>] [-r [[<strand>]<chr>:]<start>..<end> [-R]]
[-CTVNJMKQAFGUBHZWTOLE] [-w <exons.fa>] [-x <cds.fa>] [-y <tr_cds.fa>]
[-i <maxintron>]
Filters and/or converts GFF3/GTF2 records.
<input_gff> is a GFF file, use '-' if the GFF records will be given at stdin
常用参数介绍:
-g 序列文件,即GFF/GTF文件第一列ID对应的序列文件。
-i 丢弃掉内含子大于的转录本(mRNA/transcript)
-r 起始和终止位置,填写示例100.10000即为输出与100到10000有重叠的所有转录组,也可以限制序列ID及链,填写示例:+Chr1:100..10000。
-R 丢弃掉此范围的转录本,与-r相反。
-U 丢弃掉 single-exon的转录本
-C 丢低调无CDS的转录本。
-V 丢弃掉含有移码突变的转录本。
-H 如果使用了-V,则重新检查并调整内含子相位,避免由于翻译起始位点选择的位置不对导致移码突变的产生。
-B 如果使用了-V, 对于单外显子基因,则重新检查相反的链,是否存在移码突变。
-N 丢弃掉多外显子基因剪接位点不是常见的 GT-AG, GC-AG or AT-AC序列。
-J 丢弃掉没有起始密码子或者终止密码子的转录本,仅保留有完整编码框的转录本。
--no-pseudo: 过滤掉含有 'pseudo' 的注释信息
-M/--merge : 合并完全相同的或者存在包含关系的转录本。
-d:使用 -M 将合并信息输出到文件中
--cluster-only: 类似于 --merge 但是不合并转录本
-K 对于-M 选项:also collapse shorter, fully contained transcripts
with fewer introns than the container
-Q 对于-M 选项:移除包含关系的转录本的限制条件:多外显子转录本将会合并,如果他们内含子位置完全一样,单外显子转录本只需要有80%一样即可合并。
--force-exons: 使GFF features的最小层级为exon
-E 对于重复的 ID或者 GFF/GTF 其他潜在的格式问题给出警告信息。
-Z 将内含子小于4 bp的邻近的两个外显子合并为一个。
-w 输出每个转录本的外显子序列
-x 输出CDS序列
-W 对于 -w 和 -x 选项,输出外显子位置坐标到FASTA序列的ID中
-y 输出蛋白质序列
-L 将Ensembl GTF 转换为 GFF3 conversion (implies -F; should be used with -m)
-o 输出"filtered" 后的GFF文件 。
-T -o 参数将输出 GTF格式。
示例:
GFF转换GTF
gffread input.gff3 -T -o out.gtf
GTF转换GFF3
gffread input.gtf -o out.gff3
根据GFF或者GTF提取CDS,蛋白质和外显子序列
gffread gene.gff3 -g genome.fa -x cds.fa -y pep.fa -w cdna.fa
只提取翻译后蛋白序列
gffread gencode.v19.annotation.gff3 -g hg19.fa -y tr_pep.fa
根据reference提取CDS序列
gffread gencode.v19.annotation.gff3 -g hg19.fa -x cds.fa
只提取外显子序列
gffread gencode.v19.annotation.gff3 -g hg19.fa -w exons.fa
你不知道的事:
1)如果genome.fa中有N、或非ATCG字符,会在转成CDS的时候变为H。:
提取后的cds file:
2)gff 注意正负链问题:
genome:XXXXXXXXXXXXTaGXXXXXXXXXXXX
cds: XXXXXXXXXXXXXCtAXXXXXXXXXXXX
当genome中负链上的碱基有改变:某个A改变成a,这时候提取CDS序列后,会发现是T变t。
那么后面插入一个T就相当于前面插入一个a;
3)gffread是提取序列,不是生成bed文件。