ゲノムの領域を Binning したデータを作ってあれこれ
必要なもの
bedtools
hg38.info
染色体ごとの長さ(bp including gaps) を記録したテキストファイル
これ https://genome.ucsc.edu/cgi-bin/hgTracks?chromInfoPage=
染色体番号でソートしておく
たとえば 50000 ごとに切る
code:makewindows
$ bedtools makewindows -g hg38.sorted.info -w 50000 > hg38.50000.bin
そのままだと隣り合う bin 端が被ってしまうので awk で2カラム目に 1 加えるという地味なことをする
こうなる
code:binning/chr1.hg38.50000.bin
$ head binning/chr1.hg38.50000.bin
chr1 1 50000
chr1 50001 100000
chr1 100001 150000
chr1 150001 200000
chr1 200001 250000
chr1 250001 300000
chr1 300001 350000
chr1 350001 400000
chr1 400001 450000
chr1 450001 500000
これに対して MACS2 で call した ChIP-Seq の peak の情報を map したりできる
code:map
$ bedtools map -a hg38.50000.bin -b ENCFF165VDC.sorted.bed -o count
Bin x TSS
各 bin から各遺伝子のTSSまでの距離を値に持つ行列を作る
縦: bin
横: 遺伝子
セルの値: bin から遺伝子のtssまでの距離
TSSの情報は UCSC の refflat.txt.gz から取ってくる
http://hgdownload.cse.ucsc.edu/goldenpath/hg38/database/
google "ucsc hg38 refflat"
refflat には transcript の ID (RefSeq ID) が載っているのでこれを NCBIGene ID に変換する
TogoID を使って一発
重複は除去しておく
code:join
$ head refseq.tss.tsv
NM_000014 chr12 9067707 9067707 . -
NM_000015 chr8 18391281 18391281 . +
NM_000016 chr1 75724708 75724708 . +
NM_000017 chr12 120725825 120725825 . +
NM_000018 chr17 7219937 7219937 . +
NM_000019 chr11 108121566 108121566 . +
NM_000020 chr12 51907503 51907503 . +
NM_000021 chr14 73136506 73136506 . +
NM_000022 chr20 44619521 44619521 . -
NM_000023 chr17 50166004 50166004 . +
$ head refseq_rna-ncbigene.sorted.tsv
NM_000014 2
NM_000015 10
NM_000016 34
NM_000017 35
NM_000018 37
NM_000019 38
NM_000020 94
NM_000021 5663
NM_000022 100
NM_000023 6442
$ join refseq.tss.tsv refseq_rna-ncbigene.sorted.tsv | awk 'BEGIN{OFS="\t"} { print $2, $3, $4, $7, $5, $6 }' | sort -u > genes.tss.bed
$ head genes.tss.bed
chr1 100038094 100038094 64645 . +
chr1 100083569 100083569 163786 . -
chr1 100133149 100133149 54482 . +
chr1 100148447 100148447 127495 . -
chr1 100186918 100186918 1629 . -
chr1 100264741 100264741 100506007 . -
chr1 100266157 100266157 8634 . +
chr1 100266215 100266215 8634 . +
chr1 100281240 100281240 693138 . +
chr1 10032957 10032957 10277 . +
bedtools closest -d で距離が出せる
全 bin x 全遺伝子TSSまでの距離が欲しいので iterate する
遺伝子を染色体ごとに分ける
code:iterate
$ cat genes.tss.bed | cut -f 1 | sort -u | while read chr; do awk -v chr=${chr} '$1 == chr' genes.tss.bed | bedtools sort > ${chr}.genes.tss.bed; done
$ ls
chr1.genes.tss.bed chr12.genes.tss.bed chr15.genes.tss.bed chr18.genes.tss.bed chr20.genes.tss.bed chr3.genes.tss.bed chr6.genes.tss.bed chr9.genes.tss.bed chrY.genes.tss.bed
chr10.genes.tss.bed chr13.genes.tss.bed chr16.genes.tss.bed chr19.genes.tss.bed chr21.genes.tss.bed chr4.genes.tss.bed chr7.genes.tss.bed chrM.genes.tss.bed genes.tss.bed
chr11.genes.tss.bed chr14.genes.tss.bed chr17.genes.tss.bed chr2.genes.tss.bed chr22.genes.tss.bed chr5.genes.tss.bed chr8.genes.tss.bed chrX.genes.tss.bed
$ cd ../..
$ ls binning
chr1.hg38.5000.bin chr12.hg38.50000.bin chr16.hg38.5000.bin chr19.hg38.50000.bin chr22.hg38.5000.bin chr5.hg38.50000.bin chr9.hg38.5000.bin chrY.hg38.50000.bin
chr1.hg38.50000.bin chr13.hg38.5000.bin chr16.hg38.50000.bin chr2.hg38.5000.bin chr22.hg38.50000.bin chr6.hg38.5000.bin chr9.hg38.50000.bin hg38.5000.bin
chr10.hg38.5000.bin chr13.hg38.50000.bin chr17.hg38.5000.bin chr2.hg38.50000.bin chr3.hg38.5000.bin chr6.hg38.50000.bin chrM.hg38.5000.bin hg38.50000.bin
chr10.hg38.50000.bin chr14.hg38.5000.bin chr17.hg38.50000.bin chr20.hg38.5000.bin chr3.hg38.50000.bin chr7.hg38.5000.bin chrM.hg38.50000.bin
chr11.hg38.5000.bin chr14.hg38.50000.bin chr18.hg38.5000.bin chr20.hg38.50000.bin chr4.hg38.5000.bin chr7.hg38.50000.bin chrX.hg38.5000.bin
chr11.hg38.50000.bin chr15.hg38.5000.bin chr18.hg38.50000.bin chr21.hg38.5000.bin chr4.hg38.50000.bin chr8.hg38.5000.bin chrX.hg38.50000.bin
chr12.hg38.5000.bin chr15.hg38.50000.bin chr19.hg38.5000.bin chr21.hg38.50000.bin chr5.hg38.5000.bin chr8.hg38.50000.bin chrY.hg38.5000.bin
$ ls ucsc_tss/genes_tss
chr1.gene.tss.bed chr12.gene.tss.bed chr15.gene.tss.bed chr18.gene.tss.bed chr20.gene.tss.bed chr3.gene.tss.bed chr6.gene.tss.bed chr9.gene.tss.bed chrY.gene.tss.bed
chr10.gene.tss.bed chr13.gene.tss.bed chr16.gene.tss.bed chr19.gene.tss.bed chr21.gene.tss.bed chr4.gene.tss.bed chr7.gene.tss.bed chrM.gene.tss.bed genes.tss.bed
chr11.gene.tss.bed chr14.gene.tss.bed chr17.gene.tss.bed chr2.gene.tss.bed chr22.gene.tss.bed chr5.gene.tss.bed chr8.gene.tss.bed chrX.gene.tss.bed
# closest -d のテスト
$ bedtools closest -d -a binning/chr1.hg38.50000.bin -b ucsc_tss/genes_tss/chr1.genes.tss.bed | head
chr1 1 50000 chr1 11873 11873 100287102 . + 0
chr1 1 50000 chr1 14361 14361 653635 . - 0
chr1 1 50000 chr1 17368 17368 103504738 . - 0
chr1 1 50000 chr1 17368 17368 102466751 . - 0
chr1 1 50000 chr1 17368 17368 102465910 . - 0
chr1 1 50000 chr1 17368 17368 102465909 . - 0
chr1 1 50000 chr1 30365 30365 100422831 . + 0
chr1 1 50000 chr1 30365 30365 100422834 . + 0
chr1 1 50000 chr1 30365 30365 100302278 . + 0
chr1 1 50000 chr1 30365 30365 100422919 . + 0
# closest で -b を1エントリ(=1gene)のみにした場合のテスト
$ bedtools closest -d -a binning/chr1.hg38.50000.bin -b <(head -1 ucsc_tss/genes_tss/chr1.genes.tss.bed) | head
chr1 1 50000 chr1 11873 11873 100287102 . + 0
chr1 50001 100000 chr1 11873 11873 100287102 . + 38128
chr1 100001 150000 chr1 11873 11873 100287102 . + 88128
chr1 150001 200000 chr1 11873 11873 100287102 . + 138128
chr1 200001 250000 chr1 11873 11873 100287102 . + 188128
chr1 250001 300000 chr1 11873 11873 100287102 . + 238128
chr1 300001 350000 chr1 11873 11873 100287102 . + 288128
chr1 350001 400000 chr1 11873 11873 100287102 . + 338128
chr1 400001 450000 chr1 11873 11873 100287102 . + 388128
chr1 450001 500000 chr1 11873 11873 100287102 . + 438128
# 整形する
# chromosome, bin start, bin end, NCBI Gene ID, distance, strand, TSS position
$ bedtools closest -d -a binning/chr1.hg38.50000.bin -b <(head -1 ucsc_tss/genes_tss/chr1.genes.tss.bed) | awk 'BEGIN{FS=OFS="\t"} { print $1, $2, $3, $7, $10, $9, $5 }' | head
chr1 1 50000 100287102 0 + 11873
chr1 50001 100000 100287102 38128 + 11873
chr1 100001 150000 100287102 88128 + 11873
chr1 150001 200000 100287102 138128 + 11873
chr1 200001 250000 100287102 188128 + 11873
chr1 250001 300000 100287102 238128 + 11873
chr1 300001 350000 100287102 288128 + 11873
chr1 350001 400000 100287102 338128 + 11873
chr1 400001 450000 100287102 388128 + 11873
chr1 450001 500000 100287102 438128 + 11873
( bedtools closest -d -a binning/chr1.hg38.50000.bin -b | awk ; ) 0.02s user 0.08s system 66% cpu 0.161 total
# これを全bin x 全遺伝子 per chromosome する
$ for binw in 5000 50000; do cut -f 1 binning/hg38.50000.bin | sort -u | while read chr; do cat ucsc_tss/genes_tss/${chr}.genes.tss.bed | while read bedline; do bedtools closest -d -a binning/${chr}.hg38.${binw}.bin -b <(echo $bedline) | awk 'BEGIN{FS=OFS="\t"} { print $1, $2, $3, $7, $10, $9, $5 }'; done; done > hg38.${binw}.bin-gene-tss-distance.bed; done
# 当然終わらない
bedtools closest で -k オプションをつけるとmatchするbinを増やせるので全部を対象にすれば遺伝子側を分割しなくてもいいのでは？
code:closest-k
$ time (bedtools closest -k $(wc -l ucsc_tss/genes_tss/chr1.genes.tss.bed | awk '$0=$1') -d -a binning/chr1.hg38.50000.bin -b ucsc_tss/genes_tss/chr1.genes.tss.bed | wc -l )
19969800
( bedtools closest -k -d -a binning/chr1.hg38.50000.bin -b | wc -l; ) 31.38s user 1.19s system 104% cpu 31.251 total
$ time (bedtools closest -k $(wc -l ucsc_tss/genes_tss/chr1.genes.tss.bed | awk '$0=$1') -d -a binning/chr1.hg38.5000.bin -b ucsc_tss/genes_tss/chr1.genes.tss.bed | wc -l )
199665920
( bedtools closest -k -d -a binning/chr1.hg38.5000.bin -b | wc -l; ) 314.86s user 10.90s system 104% cpu 5:12.11 total
さくっと grid engine でバラす
code:gene_tss_distance.sh
#!/bin/bash
set -x
project_dir="${HOME}/repos/machima/data"
for binw in 5000 50000; do
cut -f 1 ${project_dir}/binning/hg38.50000.bin | sort -u | while read chr; do
qsub -cwd -j y -o "${project_dir}/out" ${project_dir}/script/gene_tss_distance.job.sh "${binw}" "${chr}"
done
done
code: gene_tss_distance.job.sh
#!/bin/bash
#$ -S /bin/bash
#$ -j y
binw=${1}
chr=${2}
data_dir="${HOME}/repos/machima/data"
outdir="${data_dir}/out/${binw}/${chr}"
out="${outdir}/out.bed"
mkdir -p ${outdir}
binning="${data_dir}/binning/${chr}.hg38.${binw}.bin"
genes_tss="${data_dir}/ucsc_tss/genes_tss/${chr}.genes.tss.bed"
${HOME}/local/bin/bedtools closest -k $(wc -l ${genes_tss} | awk '$0=$1') -d -a ${binning} -b ${genes_tss} |\
awk 'BEGIN{FS=OFS="\t"} { print $1, $2, $3, $7, $10, $9, $5 }' > ${out}
ちゃんと全部取れている？
code:wc
$ wc -l binning/chr1.hg38.50000.bin
4980 binning/chr1.hg38.50000.bin
$ wc -l ucsc_tss/genes_tss/chr1.genes.tss.bed
4010 ucsc_tss/genes_tss/chr1.genes.tss.bed
$ wc -l out/50000/chr1/out.bed
19969800 out/50000/chr1/out.bed
$ echo "4980 * 4010" | bc
19969800
大丈夫そう
できあがった bed ファイルのカラムは
chromosome, bin start, bin end, NCBI Gene ID, distance from bin to TSS, strand, TSS position
問題
Transcript variant がある遺伝子は TSS を複数持つ
code:chr1/out.bed
chr1 1 50000 148398 873922 + 923922
chr1 1 50000 148398 875730 + 925730
https://www.ncbi.nlm.nih.gov/gene/148398
Transcript variant 1/2 と 3 で 左端の Exon のポジションが違うので TSS も違う
とりあえずデータには全部突っ込んでいるが、、
Bin x Enhancer
各 Bin に Enhancer が含まれているか否かを 0/1 で表現する
全ての Enhancer に対してこれをやる
Bin x Gene = enhancer_included? の行列
bedtools map で overwrap が出せるので、binning に対して enhancer を map すれば OK
Enhancer の情報はSlidebase (FANTOMの系列) からダウンロードできるものを使う
https://slidebase.binf.ku.dk/human_enhancers/presets
section 1. Extensive enhancers の enhancer_tss_associations.bed をダウンロード
code:enhancer_tss_associations.bed
$ head enhancer_tss_associations.bed
chr1 858252 861621 chr1:858256-858648;NM_152486;SAMD11;R:0.404;FDR:0 404 . 858452 858453 0,0,0 2 401,1001 0,2368
chr1 894178 956888 chr1:956563-956812;NM_015658;NOC2L;R:0.202;FDR:8.01154668254404e-08 202 . 956687 956688 0,0,0 2 1001,401 0,62309
chr1 901376 956888 chr1:956563-956812;NM_001160184,NM_032129;PLEKHN1;R:0.422;FDR:0 422 . 956687 956688 0,0,0 2 1001,401 0,55111
chr1 901376 1173762 chr1:1173386-1173736;NM_001160184,NM_032129;PLEKHN1;R:0.311;FDR:0 311 . 1173561 1173562 0,0,0 2 1001,401 0,271985
chr1 935051 942164 chr1:941791-942135;NM_001142467,NM_021170;HES4;R:0.187;FDR:6.32949888009368e-07 187 . 941963 941964 0,0,0 2 1001,401 0,6712
chr1 935051 1005621 chr1:1005293-1005547;NM_001142467,NM_021170;HES4;R:0.236;FDR:6.28221217150423e-11 236 . 1005420 1005421 0,0,0 2 1001,401 0,70169
chr1 935051 1124972 chr1:1124600-1124942;NM_001142467,NM_021170;HES4;R:0.225;FDR:6.02784045803393e-10 225 . 1124771 1124772 0,0,0 2 1001,401 0,189520
chr1 935051 1136470 chr1:1136075-1136463;NM_001142467,NM_021170;HES4;R:0.27;FDR:2.85550080524545e-14 270 . 1136269 1136270 0,0,0 2 1001,401 0,201018
chr1 935051 1137107 chr1:1136667-1137146;NM_001142467,NM_021170;HES4;R:0.164;FDR:8.84099081246056e-06 164 . 1136906 1136907 0,0,0 2 1001,401 0,201655
chr1 945701 1335218 chr1:945769-946034;NM_001039577,NM_030937;CCNL2;R:0.176;FDR:3.21786973606028e-06 176 . 945901 945902 0,0,0 2 401,1001 0,388516
必要なのは NCBI gene ID に対する enhancer region の mapping なので、まずは 4th column をパースする
code:awk
$ cat enhancer_tss_associations.bed | awk 'BEGIN{FS=OFS="\t"} { split($4, info, ";"); print $1, $2, $3, info2, $5, $6, info3 }' | head
chr1 858252 861621 NM_152486 404 . SAMD11
chr1 894178 956888 NM_015658 202 . NOC2L
chr1 901376 956888 NM_001160184,NM_032129 422 . PLEKHN1
chr1 901376 1173762 NM_001160184,NM_032129 311 . PLEKHN1
chr1 935051 942164 NM_001142467,NM_021170 187 . HES4
chr1 935051 1005621 NM_001142467,NM_021170 236 . HES4
chr1 935051 1124972 NM_001142467,NM_021170 225 . HES4
chr1 935051 1136470 NM_001142467,NM_021170 270 . HES4
chr1 935051 1137107 NM_001142467,NM_021170 164 . HES4
chr1 945701 1335218 NM_001039577,NM_030937 176 . CCNL2
ここでも複数のTranscriptに紐付いているが1行は1遺伝子に対応しているので最初の Transcript を NCBI Gene ID に変換する
code:awk
$ cat enhancer_tss_associations.bed | awk 'BEGIN{FS=OFS="\t"} { split($4, info, ";"); split(info2, tx, ","); print $1, $2, $3, tx1, $5, $6 }' > enhancer_transcript_mapping.bed
$ head enhancer_transcript_mapping.bed
chr1 858252 861621 NM_152486 404 .
chr1 894178 956888 NM_015658 202 .
chr1 901376 956888 NM_001160184 422 .
chr1 901376 1173762 NM_001160184 311 .
chr1 935051 942164 NM_001142467 187 .
chr1 935051 1005621 NM_001142467 236 .
chr1 935051 1124972 NM_001142467 225 .
chr1 935051 1136470 NM_001142467 270 .
chr1 935051 1137107 NM_001142467 164 .
chr1 945701 1335218 NM_001039577 176 .
# Join がコケるので sort -V で natural sort して NM_XXX を並べておく
$ sort -V ../ucsc_tss/refseq_rna-ncbigene.sorted.tsv > refseq_rna-ncbigene.sorted.tsv
$ cat enhancer_transcript_mapping.bed | awk 'BEGIN{FS=OFS="\t"} { print $4, $1, $2, $3, $5, $6 }' | sort -V > enhancer_transcript_mapping.sorted.txt
$ join enhancer_transcript_mapping.sorted.txt refseq_rna-ncbigene.sorted.tsv | awk 'BEGIN{OFS="\t"} { print $2, $3, $4, $7, $5, $6, $1 }' > enhancer_gene_mapping.bed
$ head enhancer_gene_mapping.bed
chr8 18248254 18524115 10 215 . NM_000015
chr1 76189542 76250968 34 362 . NM_000016
chr12 120665509 121164071 35 229 . NM_000017
chr11 107671444 107992758 38 208 . NM_000019
chr12 52291022 52301702 94 595 . NM_000020
chr12 52299778 52301702 94 770 . NM_000020
chr12 52300701 52419009 94 266 . NM_000020
chr12 52300701 52588132 94 257 . NM_000020
chr14 73602642 73695547 5663 603 . NM_000021
chr14 73602642 73825217 5663 424 . NM_000021
bedtools map していく
code:bedtools
$ cd ..
$ bedtools sort -i slidebase/enhancer_gene_mapping.bed > slidebase/enhancer_gene_mapping.sorted.bed
結果を見ると 50000 bin に跨って存在する enhancer がある
そんなでかいことある？
元のbedを見てもめちゃくちゃレンジが広い
なんなんだこのbedは
https://www.biostars.org/p/417137/
案の定困っている人がいた
どうもこいつは UCSC genome browser で表示するためのものであり、最初の3カラムは表示領域 (enhancer+TSS) を示しているっぽい
4th column の genomic range が enhancer 領域を示している模様、、
やりなおし
code:awk
$ cat enhancer_tss_associations.bed | awk 'BEGIN{FS=OFS="\t"} { split($4, info, ";"); split(info1, chr, ":"); split(chr2, range, "-"); split(info2, tx, ","); print chr1, range1, range2, tx1, $5, $6 }' > enhancer_transcript_mapping.bed
$ head enhancer_transcript_mapping.bed
chr1 858256 858648 NM_152486 404 .
chr1 956563 956812 NM_015658 202 .
chr1 956563 956812 NM_001160184 422 .
chr1 1173386 1173736 NM_001160184 311 .
chr1 941791 942135 NM_001142467 187 .
chr1 1005293 1005547 NM_001142467 236 .
chr1 1124600 1124942 NM_001142467 225 .
chr1 1136075 1136463 NM_001142467 270 .
chr1 1136667 1137146 NM_001142467 164 .
chr1 945769 946034 NM_001039577 176 .
$ join <(awk 'BEGIN{FS=OFS="\t"} { print $4, $0 }' enhancer_transcript_mapping.bed| sort -V) refseq_rna-ncbigene.sorted.tsv | head
NM_000015 chr8 18523718 18524110 NM_000015 215 . 10
NM_000016 chr1 76250721 76250813 NM_000016 362 . 34
NM_000017 chr12 120665661 120665757 NM_000017 229 . 35
NM_000019 chr11 107671635 107671653 NM_000019 208 . 38
NM_000020 chr12 52291135 52291309 NM_000020 595 . 94
NM_000020 chr12 52299885 52300072 NM_000020 770 . 94
NM_000020 chr12 52418719 52418897 NM_000020 266 . 94
NM_000020 chr12 52587798 52588065 NM_000020 257 . 94
NM_000021 chr14 73695241 73695452 NM_000021 603 . 5663
NM_000021 chr14 73824884 73825148 NM_000021 424 . 5663
$ join <(awk 'BEGIN{FS=OFS="\t"} { print $4, $0 }' enhancer_transcript_mapping.bed| sort -V) refseq_rna-ncbigene.sorted.tsv > enhancer_gene_mapping.bed
code:diff
$ wc -l enhancer_gene_mapping.bed
64842 enhancer_gene_mapping.bed
$ wc -l enhancer_transcript_mapping.bed
66943 enhancer_transcript_mapping.bed
$ echo 64842-66943 | bc
-2101
$ grep NM_207647 refseq_rna-ncbigene.sorted.tsv
ない
試しにdiffで出てくるものを見てみると withdrawn record だった
https://www.ncbi.nlm.nih.gov/nuccore/NM_207647
なので無視する
やっと再マッピングできます
code:bedtools
$ bedtools sort -i slidebase/enhancer_gene_mapping.bed > slidebase/enhancer_gene_mapping.sorted.bed
$ bedtools map -o distinct -c 4 -a binning/chr1.hg38.50000.bin -b slidebase/enhancer_gene_mapping.sorted.bed
今度は見た感じ違和感がない
ではこれを 5000, 50000 bins の各chromosomeに対してやります
code:for
$ for binw in 5000 50000; do cut -f 1 slidebase/enhancer_gene_mapping.sorted.bed|sort -u|while read chr; do bedtools map -o distinct -c 4 -a binning/${chr}.hg38.${binw}.bin -b slidebase/enhancer_gene_mapping.sorted.bed > ${binw}.map_enhancer.${chr}.bed; done; done
$ head 50000.map_enhancer.chr10.bed
chr10 1 50000 .
chr10 50001 100000 .
chr10 100001 150000 .
chr10 150001 200000 .
chr10 200001 250000 .
chr10 250001 300000 .
chr10 300001 350000 .
chr10 350001 400000 .
chr10 400001 450000 10771,22982,347688
chr10 450001 500000 .
できました
結論
bedtoolsはマジで偉大
ファイルについて
code:shell
$ ls data
binning refflat script slidebase tss_distance ucsc_tss
binning
bedtools makewindows の出力を染色体ごとに分けた bin データ。bin width 5000 と 50000 がある。
refflat
UCSC から落としてきたアノテーションファイル。TSSの情報に使う。
script
bin to TSS distance の計算に使ったスクリプト。遺伝研スパコンでしか動かないので可搬性ガン無視でハードコードしている。
slidebase
bin to enhancer の計算元のデータと出力データ。 slidebase/out に binwidth, chromosome ごとに分割されたデータが入っている。
tss_distance
tss_distance/5000, tss_distance/50000 それぞれに染色体ごとの bin to TSS distance のデータが入っている。ファイル名もうちょっとちゃんとすればよかった。
ucsc_tss
refflat から作った TSS の情報の bedfile が入っている。tss_distance slidebase の計算に使う。TogoID から取ってきた RefSeq to NCBI Gene の ID conversion table も入っている。
使うデータ
Bin to TSS distance
tss_distance/5000/chr1/out.bed のように binwidth, chromosome ごとのデータが入っている
カラム
1. Chromosome
2. Bin start
3. Bin end
4. NCBI Gene ID
e.g. https://identifiers.org/ncbigene/100287102
5. Distance to TSS (bp)
6. Strand
このデータではあまり意味がない
7. TSS location
bin の数 x Transcript なのでデータが無駄にバカでかい
chr1 5000 binw で 49792, 50000 binw で 4980
chr 1 の TSSの数が 4010
5000 binw で 199,665,920, 50000 binw で 19,969,800 とかになる (#lines)
ファイルサイズの合計が tss_distance 以下だけで 60GB くらいある
code:out.bed
$ head tss_distance/5000/chr1/out.bed
chr1 1 5000 100287102 6873 + 11873
chr1 1 5000 653635 9361 - 14361
chr1 1 5000 103504738 12368 - 17368
chr1 1 5000 102466751 12368 - 17368
chr1 1 5000 102465910 12368 - 17368
chr1 1 5000 102465909 12368 - 17368
chr1 1 5000 100422831 25365 + 30365
chr1 1 5000 100422834 25365 + 30365
chr1 1 5000 100302278 25365 + 30365
chr1 1 5000 100422919 25365 + 30365
Bin to Enhancer
slidebase/out/50000.map_enhancer.chr10.bed のように binwidth, chr に従ってファイル名がついている
カラムは Chromosome, Bin start, Bin end に続いて4カラム目に「そのbinの中にenhancerを持つ遺伝子のNCBI Gene ID」がcsvで入っている
code:bed
$ head slidebase/out/50000.map_enhancer.chr10.bed
chr10 1 50000 .
chr10 50001 100000 .
chr10 100001 150000 .
chr10 150001 200000 .
chr10 200001 250000 .
chr10 250001 300000 .
chr10 300001 350000 .
chr10 350001 400000 .
chr10 400001 450000 10771,22982,347688
chr10 450001 500000 .