gtfファイルの取り扱い

#gftファイルの取り扱いについて
【目的】RNA-seqの解析をしていて、featureCountsでの集計が、ensembl idになってしまって、Gene symbolなどがわからなくなってしまうので、追加でデータを挿れたい。
(他にいい方法がある場合、教えて下さい…)
(pandasの勉強も兼ねてやっているので、遠回りしているところもあります。。。)

##*gtfファイル*
今回は、GENCODEのhuman Comprehensive gene annotation のAllを使用している。
ここからダウンロード
https://www.gencodegenes.org/human/

ダウンロードすると、
gencode.v34.chr_patch_hapl_scaff.annotation.gtf
が手に入る。

##*gtfparse*
gtfの取り扱いについて調べたところ、これがヒットしたので、取り敢えず、インストールしてみた。

```Bash
pip install gtfparse
```

これでデータを見ることが出来るか確かめてみる。

```Python
import pandas as pd
from gtfparse import read_gtf

df = read_gtf("gencode.v34.chr_patch_hapl_scaff.annotation.gtf")
print(df.head())
print(df.shape)
print(df.columns)

```

これがやたら時間がかかる・・・

```Bash
INFO:root:Extracted GTF attributes: ['gene_id', 'gene_type', 'gene_name', 'level', 'hgnc_id', 'havana_gene', 'transcript_id', 'transcript_type', 'transcript_name', 'transcript_support_level', 'tag', 'havana_transcript', 'exon_number', 'exon_id', 'ont', 'protein_id', 'ccdsid']
  seqname  source     feature  start    end  score  ...     havana_transcript exon_number            exon_id ont protein_id ccdsid
0    chr1  HAVANA        gene  11869  14409    NaN  ...                                                                           
1    chr1  HAVANA  transcript  11869  14409    NaN  ...  OTTHUMT00000362751.1                                                     
2    chr1  HAVANA        exon  11869  12227    NaN  ...  OTTHUMT00000362751.1           1  ENSE00002234944.1                      
3    chr1  HAVANA        exon  12613  12721    NaN  ...  OTTHUMT00000362751.1           2  ENSE00003582793.1                      
4    chr1  HAVANA        exon  13221  14409    NaN  ...  OTTHUMT00000362751.1           3  ENSE00002312635.1                      

[5 rows x 25 columns]
(3171209, 25)
Index(['seqname', 'source', 'feature', 'start', 'end', 'score', 'strand',
       'frame', 'gene_id', 'gene_type', 'gene_name', 'level', 'hgnc_id',
       'havana_gene', 'transcript_id', 'transcript_type', 'transcript_name',
       'transcript_support_level', 'tag', 'havana_transcript', 'exon_number',
       'exon_id', 'ont', 'protein_id', 'ccdsid'],
      dtype='object')

```

こんな感じでデータが取れたので、OKそう。31万以上のデータが入ってるんだね。(Gene symbolとしては同一のものがたくさん含まれてます。)

ここで必要そうなものは、
'gene_name', 'seqname','start','end','strand','gene_type', 'gene_id','transcript_id'
なので、これらを取ってきたら、まずはOK。あたらにcsvファイルとして保存しておく。

```Python:Gene_Symbol_from_GTF.py
from gtfparse import read_gtf

df = read_gtf("gencode.v34.chr_patch_hapl_scaff.annotation.gtf")
print(df.head())
print(df.shape)
print(df.columns)

df2 = df[['gene_name', 'seqname','start','end','strand','gene_type', 'gene_id','transcript_id']]
df2.to_csv("gencode.v34.chr_patch_hapl_scaff.annotation.cur.gtf.csv")
```

一応確認。


```Bash
INFO:root:Extracted GTF attributes: ['gene_id', 'gene_type', 'gene_name', 'level', 'hgnc_id', 'havana_gene', 'transcript_id', 'transcript_type', 'transcript_name', 'transcript_support_level', 'tag', 'havana_transcript', 'exon_number', 'exon_id', 'ont', 'protein_id', 'ccdsid']
  gene_name seqname  start    end strand                           gene_type            gene_id      transcript_id
0   DDX11L1    chr1  11869  14409      +  transcribed_unprocessed_pseudogene  ENSG00000223972.5                   
1   DDX11L1    chr1  11869  14409      +  transcribed_unprocessed_pseudogene  ENSG00000223972.5  ENST00000456328.2
2   DDX11L1    chr1  11869  12227      +  transcribed_unprocessed_pseudogene  ENSG00000223972.5  ENST00000456328.2
3   DDX11L1    chr1  12613  12721      +  transcribed_unprocessed_pseudogene  ENSG00000223972.5  ENST00000456328.2
4   DDX11L1    chr1  13221  14409      +  transcribed_unprocessed_pseudogene  ENSG00000223972.5  ENST00000456328.2
5   DDX11L1    chr1  12010  13670      +  transcribed_unprocessed_pseudogene  ENSG00000223972.5  ENST00000450305.2
6   DDX11L1    chr1  12010  12057      +  transcribed_unprocessed_pseudogene  ENSG00000223972.5  ENST00000450305.2
7   DDX11L1    chr1  12179  12227      +  transcribed_unprocessed_pseudogene  ENSG00000223972.5  ENST00000450305.2
8   DDX11L1    chr1  12613  12697      +  transcribed_unprocessed_pseudogene  ENSG00000223972.5  ENST00000450305.2
9   DDX11L1    chr1  12975  13052      +  transcribed_unprocessed_pseudogene  ENSG00000223972.5  ENST00000450305.2
          gene_name     seqname   start     end strand       gene_type            gene_id      transcript_id
3171199  AC007325.2  KI270734.1  138480  138482      -  protein_coding  ENSG00000277196.4  ENST00000621424.4
3171200  AC007325.2  KI270734.1  161689  161852      -  protein_coding  ENSG00000277196.4  ENST00000621424.4
3171201  AC007325.2  KI270734.1  161587  161626      -  protein_coding  ENSG00000277196.4  ENST00000621424.4
3171202  AC007325.2  KI270734.1  138082  138482      -  protein_coding  ENSG00000277196.4  ENST00000621424.4
3171203  AC237676.1  KI270744.1   51009   51114      -           snRNA  ENSG00000278625.1                   
3171204  AC237676.1  KI270744.1   51009   51114      -           snRNA  ENSG00000278625.1  ENST00000616830.1
3171205  AC237676.1  KI270744.1   51009   51114      -           snRNA  ENSG00000278625.1  ENST00000616830.1
3171206  AC116622.1  KI270750.1  148668  148843      +           snRNA  ENSG00000277374.1                   
3171207  AC116622.1  KI270750.1  148668  148843      +           snRNA  ENSG00000277374.1  ENST00000612925.1
3171208  AC116622.1  KI270750.1  148668  148843      +           snRNA  ENSG00000277374.1  ENST00000612925.1
```

OKそうですね。(transcript_idは今回使わないのですが、ちょっと別のことで使用したいので、取っておきました。）

今回は、gene_nameとgene_type, gene_idの対応表がほしいので、gene_nameで重複の削除をします。

```Python
import pandas as pd

df = pd.read_csv('gencode.v34.chr_patch_hapl_scaff.annotation.cur.gtf.csv',header = 0, sep=',')
df2 = df[~df.duplicated(['gene_name'])]
df2.to_csv("gencode.v34.Symbol.csv")
```

これで、目的のものが持ってこれました。
あとは、自分の集計データとこれをExcelで言うところのvlookupをしたらいいです。


＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿
もしかして、そもそもgtfファイルはread_csvで読み込めた？