hg38のfastaとgtf/gff3 - mecobalamin’s diary

GRCh38とhg38はヒトゲノムのシークエンスデータ
ほとんど一緒らしい。
前にも引用させてもらったけどここに説明がある。
GRCh38とhg38の違い(含むミトコンドリア)

GRCh38染色体名(ヘッダ行)が長く(gi|568815597|ref|NC_000001.11|)、hg38は簡潔(chr1)。

ヘッダ行に違いがあるのはやや問題で、
JBrowseで表示させたりfeatureCountsでリードカウントするときに
スムーズにできないことがあった

そこで簡潔な表示で揃えてfastaとgtf/gff3を作成する

まずはシークエンスデータのダウンロード
UCSC Genome Browser Downloads

リンクをたどるといくつかファイルが見つかる
chromosomeごとに分かれているデータを使う
Index of /goldenPath/hg38/bigZips

ファイルの説明はこちら

hg38.chromFa.tar.gz - The assembly sequence in one file per chromosome.
Repeats from RepeatMasker and Tandem Repeats Finder (with period
of 12 or less) are shown in lower case; non-repeating sequence is
shown in upper case.

ダウンロードと展開

$ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.chromFa.tar.gz
$ gzip -d hg38.chromFa.tar.gz
$ tar -xvf hg38.chromFa.tar

展開するとパッチを含めて450ファイルぐらいある
パッチは取り扱いがよくわからないので
とりあえずchromosomeのファイルだけマージして使う

$ :> merge_chrom.fa
$ for i in {1..22}; do cat chr${i}.fa >> merge_chrom.fa; done
$ cat chrX.fa >> merge_chrom.fa
$ cat chrY.fa >> merge_chrom.fa

このfastaファイルをhisat2-build及びhisat2に使う
mecobalamin.hatenablog.com

次にhg38のgtfファイルを作成する
Gene annotation データを用意する（gtf形式） - Palmsonntagmorgen

NCBIからダウンロードできるgffファイルは詳しい表記のヘッダなので、
UCSCのサイトからgtfファイルをダウンロードしてgff3に変換する

Table Browser@UCSC
Table Browser

設定はこうした・・・ハズ。
clade: Mammal
genome: Human
assembry: Dec. 2013 (GRch38/hg38)
group: Gene and Gene Predictions
track: NCBI RefSeq
table: RefSeq All (ncbiRefSeq)
region: genome
output format: GTF - gene transfer format (limited)
output file: UCSC.hg38.gtf.gz
file type returned: gzip compressed

gtfファイルをgffreadをつかってgff3に変換する

$ gffread -E UCSC.hg38.gtf -o- > UCSC.hg38.gff3

このファイルはパッチで修正？された配列を含む
つまりパッチとchromosomeの間に重複がある
正しいやり方かよくわからないがパッチ側の重複した配列を除く
取り除かないと通らないコマンドがあった
bam2wigだったかな

$ cat UCSC.hg38.gff3 | awk '(a[$9]++ < 1000) && ($1 !~/[_alt]/) {print}' > rename_UCSC.hg38.gff3

exonやCDSで同じ名前のエントリーがあるので
1000個以上の重複で"_alt"を含む場合に除く
コードでは重複が1000個以下でかつ"_alt"を含まない行をファイルに出力している
1000個以上重複している配列はなさそうだったので、
使った感じでは取りこぼしはなさそう

今度はgff3をgtfに変換する

$  gffread rename_UCSC.hg38.gff3 -T -o- > rename_UCSC.hg38.gtf

gtfはfeatureCountsに、gff3はJBrowseに使う

ファイルの中身はこんな感じ
rename_UCSC.hg38.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";

rename_UCSC.hg38.gff3

# gffread -E UCSC.hg38.gtf -o-
# gffread v0.10.6
##gff-version 3
chr1    hg38_ncbiRefSeq transcript      11874   14409   .       +       .       ID=NR_046018.2;geneID=NR_046018.2
chr1    hg38_ncbiRefSeq exon    11874   12227   .       +       .       Parent=NR_046018.2
chr1    hg38_ncbiRefSeq exon    12613   12721   .       +       .       Parent=NR_046018.2
chr1    hg38_ncbiRefSeq exon    13221   14409   .       +       .       Parent=NR_046018.2
chr1    hg38_ncbiRefSeq transcript      14362   29370   .       -       .       ID=NR_024540.1;geneID=NR_024540.1
chr1    hg38_ncbiRefSeq exon    14362   14829   .       -       .       Parent=NR_024540.1
chr1    hg38_ncbiRefSeq exon    14970   15038   .       -       .       Parent=NR_024540.1