得到一個物種所有基因的TSS(轉錄起始位點)區域的bed檔案。
阿新 • • 發佈:2022-05-03
首先在UCSC的table browser 裡面下載下面這個檔案:
可以看到我這裡選擇的mm10的refseq系統的所有基因,共有29037個不同的tss,36872個轉錄本,只有24540個基因,說明有部分基因有多個tss,這個其實挺麻煩的。
#bin name chrom strand txStart txEnd cdsStart cdsEnd exonCount exonStarts exonEnds score name2 cdsStartStat cdsEndStat exonFrames0 NM_001282945 chr1 - 134199214 134235457 134202950 134234355 3 134199214,134234014,134235227, 134203590,134234446,134235457, 0 Adora1 cmpl cmpl 2,0,-1,0 NM_001039510 chr1 - 134199214 134235457 134202950 134234355 3 134199214,134234014,134235227, 134203590,134234412,134235457, 0 Adora1 cmpl cmpl 2,0,-1,0 NM_001291930 chr1 - 134199214 134235457 134202950 134203505 2 134199214,134235227, 134203590,134235457, 0 Adora1 cmpl cmpl 0,-1,0 NM_001291928 chr1 - 134199214 134234856 134202950 134234733 2 134199214,134234662, 134203590,134234856, 0 Adora1 cmpl cmpl 2,0,0 NM_001008533 chr1 - 134199214 134235457 134202950 134234355 2 134199214,134234014, 134203590,134235457, 0 Adora1 cmpl cmpl 2,0,
其實裡面可以設定直接下載所有基因的TSS區域的bed檔案,可是我不會設定各種引數,也懶得去摸索,直接對上面的檔案我可以寫指令碼處理得到需要的資料形式。 需要輸出的是bed格式檔案,如下: chrom / chromStart /chromEnd /name /score /strand 我這裡定義的TSS(轉錄起始位點)區域上下游2.5kb,所以程式碼如下:
perl -alne '{next if /^#/;if($F[3] eq "+"){$start=$F[4]-2500;$end=$F[4]+2500}else{$start=$F[5]-2500;$end=$F[5]+2500}print join("t",$F[2],$start,$end,$F[12],0,$F[3])}' ucsc.refseq.txt |sort -u >ucsc.refseq.tss.bed
最後得到的檔案如下:
tail ucsc.refseq.tss.bed chrY 816212 821212 Uba1y 0 +chrY 81798997 81803997 Gm20747 0 -chrY 82222714 82227714 Gm20736 0 +chrY 83925411 83930411 Gm20854 0 -chrY 85527019 85532019 Gm20854 0 -chrY 8832669 8837669 Gm20815 0 -chrY 895287 900287 Kdm5d 0 +chrY 90752550 90757550 G530011O06Rik 0 -chrY 90782941 90787941 Erdr1 0 +chrY 90836906 90841906 G530011O06Rik 0
這裡面會有一個問題,對於部分基因在非正常染色體的,會出現如下詭異的結果,建議乾脆刪除掉。
chr4_GL456216_random 13380 18380 Dhrsx 0 +chr4_GL456350_random -1369 3631 Ccl21c 0 -chr4_GL456350_random -1369 3631 Gm10591 0 -chr4_GL456350_random -1369 3631 Gm13304 0 -
記住,這個時候,部分基因還有多個tss哦,反正取決於你的下游分析流程啦。