フォーマット変換、GTF -> BED -> BAM

GTFファイルからBAMファイルにフォーマットを変換したい

まずフォーマットについて。
GTF/GFF
GTFとGFFフォーマット - macでインフォマティクス

SAM/BAM
BAM – 遺伝子発現解析（マイクロアレイ解析, RNA-seq）
https://bi.biopapyrus.jp/rnaseq/mapping/sam.html

あとgoogleの検索で見つかったpdfファイルが詳しかった
"NGS, 基本データフォーマット"で検索かけると
基礎生物学研究所での講義ファイルが見つかる
PDFファイルなのでリンクは貼らない

次にコマンドについて

BEDOPSはGTF/GFFをBEDに変換できる
BEDOPS: the fast, highly scalable and easily-parallelizable genome analysis toolkit — BEDOPS v2.4.41
VCF, GTF, GFF などを BED に変換する BEDOPS - macでインフォマティクス

インストールはlinux用バイナリをダウンロードして展開
パスの通っているディレクトリにリンクを貼る
使っている環境ではバイナリを~/usr/localにおいて
リンクを~/binに張っている

convert2bedとsort-bedはgtf2bedの実行に必要

$ ln -s $HOME/usr/local/bedops/bin/gff2bed $HOME/bin
$ ln -s $HOME/usr/local/bedops/bin/gtf2bed $HOME/bin
$ ln -s $HOME/usr/local/bedops/bin/convert2bed $HOME/bin
$ ln -s $HOME/usr/local/bedops/bin/sort-bed $HOME/bin

BEDからBAMへはbedtoolsでできる
https://cell-innovation.nig.ac.jp/surfers/BEDtools.html
https://bi.biopapyrus.jp/rnaseq/mapping/format-convert.html
ここからダウンロード
Installation — bedtools 2.30.0 documentation
展開してパスの通ったディレクトリにリンクを貼る

準備できたのでさっそく使ってみる

まずgtf2bedを使ってhg38のgtfをbedファイルに変換する
次にbedtoossを使ってbedをbamに変換する
bedtoolsのディレクトリにあるhuman.hg38.genomeをgenomeサイズとして指定する

$ gtf2bed < rename_UCSC.hg38.gtf > rename_UCSC.hg38.bed
$ bedtools bedtobam -i rename_UCSC.hg38.bed -g $HOME/usr/local/bedtools2/genomes/human.hg38.genome > rename_UCSC.hg38.bam

bashのリダイレクトの使い方未だによくわかんないんだけど
gtfファイルをgtf2bedに渡してbedファイルに出力しているらしい

変換前のgtfファイルはこうで

chr1    hg38_ncbiRefSeq exon    11874   12227   .       +       .       transcript_id "NR_046018.2"; gene_id "NR_046018.2";
chr1    hg38_ncbiRefSeq exon    12613   12721   .       +       .       transcript_id "NR_046018.2"; gene_id "NR_046018.2";
chr1    hg38_ncbiRefSeq exon    13221   14409   .       +       .       transcript_id "NR_046018.2"; gene_id "NR_046018.2";
chr1    hg38_ncbiRefSeq exon    14362   14829   .       -       .       transcript_id "NR_024540.1"; gene_id "NR_024540.1";
chr1    hg38_ncbiRefSeq exon    14970   15038   .       -       .       transcript_id "NR_024540.1"; gene_id "NR_024540.1";
chr1    hg38_ncbiRefSeq exon    15796   15947   .       -       .       transcript_id "NR_024540.1"; gene_id "NR_024540.1";

変換後のBAMファイルはこう

NR_046018.2     0       chr1    11874   255     354M    *       0       0       *       *
NR_046018.2     0       chr1    12613   255     109M    *       0       0       *       *
NR_046018.2     0       chr1    13221   255     1189M   *       0       0       *       *
NR_024540.1     16      chr1    14362   255     468M    *       0       0       *       *
NR_024540.1     16      chr1    14970   255     69M     *       0       0       *       *
NR_024540.1     16      chr1    15796   255     152M    *       0       0       *       *