2023.2.24 追記
1. bedtoolsによるプロモーター配列の抽出
言わずと知れた「macでインフィマティクス」(https://kazumaxneo.hatenablog.com/entry/2017/07/24/221538) およびBiostarsの記事 (https://www.biostars.org/p/17162/) を参考にしつつ、より細かい手順を示した。
まずfastaのインデックス（.fai）を利用して染色体サイズを書いたファイルを準備する。今回は、Phytozomeからダウンロードしたソルガム（Sorghum bicolor）のリファレンスゲノムを用いた。
samtools faidx Sbicolor_454_v3.1.1.fa
いったんファイルをチェックしてみる。
Chr01 80884392 7 60 61
Chr02 77742459 82232480 60 61
Chr03 74386277 161270654 60 61
Chr04 68658214 236896710 60 61
Chr05 71854669 306699235 60 61
Chr06 61277060 379751489 60 61
Chr07 65505356 442049841 60 61
Chr08 62686529 508646960 60 61
Chr09 59416394 572378272 60 61
Chr10 61233695 632784947 60 61
super_16 1391575 695039214 60 61
super_18 1172180 696453992 60 61
super_20 888508 697645719 60 61
super_22 639901 698549046 60 61
１列目が染色体名、２列目が各染色体のサイズ（bp）なので、cutでこの２列を抽出する。なお、11行目以降は不要なので、headで対応する。
cut -f 1,2 reference.fa.fai | head -10 > genome.txt
中身を確認。
Chr01 80884392
Chr02 77742459
Chr03 74386277
Chr04 68658214
Chr05 71854669
Chr06 61277060
Chr07 65505356
Chr08 62686529
Chr09 59416394
Chr10 61233695
次に、ターゲット遺伝子のリストを作成する。スタートはgff3ファイルとした。Phytozomeなどから入手できます。
##gff-version 3
##annot-version v3.1.1
##species Sorghum bicolor
Chr01 phytozomev12 gene 1951 2616 . + . ID=Sobic.001G000100.v3.2;Name=Sobic.001G000100;ancestorIdentifier=Sobic.001G000100.v2.1
Chr01 phytozomev12 mRNA 1951 2616 . + . ID=Sobic.001G000100.1.v3.2;Name=Sobic.001G000100.1;pacid=37943753;longest=1;ancestorIdentifier=Sobic.001G000100.1.v2.1;Parent=Sobic.001G000100.v3.2
Chr01 phytozomev12 CDS 1951 2454 . + 0 ID=Sobic.001G000100.1.v3.2.CDS.1;Parent=Sobic.001G000100.1.v3.2;pacid=37943753
Chr01 phytozomev12 CDS 2473 2616 . + 0 ID=Sobic.001G000100.1.v3.2.CDS.2;Parent=Sobic.001G000100.1.v3.2;pacid=37943753
Chr01 phytozomev12 gene 11180 14899 . - . ID=Sobic.001G000200.v3.2;Name=Sobic.001G000200;ancestorIdentifier=Sobic.001G000200.v2.1
Chr01 phytozomev12 mRNA 11180 14899 . - . ID=Sobic.001G000200.1.v3.2;Name=Sobic.001G000200.1;pacid=37944019;longest=1;ancestorIdentifier=Sobic.001G000200.1.v2.1;Parent=Sobic.001G000200.v3.2
Chr01 phytozomev12 CDS 14503 14601 . - 0 ID=Sobic.001G000200.1.v3.2.CDS.1;Parent=Sobic.001G000200.1.v3.2;pacid=37944019
以下のような、ターゲット遺伝子のリストも用意しておく。今回は適当な6つの遺伝子のみだが、RNA-seq解析からえられたDEGs (Differentially Expressed Genes) などを想定しています。
Sobic.005G047100
Sobic.003G270300
Sobic.001G192800
Sobic.010G101900
Sobic.002G374700
Sobic.005G035100
次に、以下のスクリプトを用いて情報を抽出・整理します。
import re
#####
my_genes_list = './my_genes.txt'
input_gff = './Sbicolor_454_v3.1.1.gene.gff3'
output_bed = './my_genes.bed'
#####
with open(my_genes_list, mode='r') as inp:
gene_list = inp.read().splitlines()
with open(input_gff, mode = 'r') as inp:
with open(output_bed, mode = 'w') as out:
for line in inp:
if not line.startswith('##'):
line = line.rstrip('\n|\r|\r\n')
line = line.split('\t')
if (line[2] == 'gene') & line[0].startswith('Chr'):
Chr = line[0]
Start = line[3]
End = line[4]
Strand = line[6]
gene_id_pattern = re.compile('Sobic\.0(0|1)[0-9]G[0-9]{6}')
gene_id = gene_id_pattern.search(line[8]).group()
if gene_id in gene_list:
out.write(f'{Chr}\t{Start}\t{End}\t{gene_id}\t0\t{Strand}\n')
else:
pass
else:
pass
else:
pass
出力されたmy_genes.bedは、こんな感じ。
Chr01 17131792 17134734 Sobic.001G192800 0 +
Chr02 73253684 73259983 Sobic.002G374700 0 +
Chr03 60683208 60688185 Sobic.003G270300 0 -
Chr05 3149924 3156356 Sobic.005G035100 0 +
Chr05 4469547 4472327 Sobic.005G047100 0 +
Chr10 9424298 9427933 Sobic.010G101900 0 -
2023.2.24 追記 ---ここから---
Twitterにて、gffreadを用いた"別解"を教えていただいた。ここでは、my_genes.txtに書く遺伝子IDをgff3内のID=に一致させる必要がある。
Sobic.005G047100.v3.2
Sobic.003G270300.v3.2
Sobic.001G192800.v3.2
Sobic.010G101900.v3.2
Sobic.002G374700.v3.2
Sobic.005G035100.v3.2
ここからgffreadとgff2bedにより、
gffread -O ./Sbicolor_454_v3.1.1.gene.gff3 --ids ./my_genes_2.txt | gff2bed > my_genes_2.bed
今回は各transcriptには興味がないため、-Oオプションでnon-transcriptのみを対象としている。このオプションを使うと--bedが使えないという"罠"があるので、gff2bedに渡してbedを作っている。
Chr01 17131791 17134734 Sobic.001G192800.v3.2 . + phytozomev12 gene . ID=Sobic.001G192800.v3.2
Chr02 73253683 73259983 Sobic.002G374700.v3.2 . + phytozomev12 gene . ID=Sobic.002G374700.v3.2
Chr03 60683207 60688185 Sobic.003G270300.v3.2 . - phytozomev12 gene . ID=Sobic.003G270300.v3.2
Chr05 3149923 3156356 Sobic.005G035100.v3.2 . + phytozomev12 gene . ID=Sobic.005G035100.v3.2
Chr05 4469546 4472327 Sobic.005G047100.v3.2 . + phytozomev12 gene . ID=Sobic.005G047100.v3.2
Chr10 9424297 9427933 Sobic.010G101900.v3.2 . - phytozomev12 gene . ID=Sobic.010G101900.v3.2
2023.2.24 追記 ---おわり---
ここから、bedtoolsで、
bedtools flank -i my_genes.bed -g genome.txt -l 2000 -r 0 -s > genes_2kb_promoters.bed
として、プロモーター領域のbedファイルを作成する。-sオプションは鎖の方向を考慮する上で必須。ないと、negative strandにコードされた遺伝子からは下流を拾ってきてしまう。最後に、refarenceのfastaファイルをもとに実際の塩基配列のmulti fasta を生成する。
bedtools getfasta -fi Sbicolor_454_v3.1.1.fa -bed genes_2kb_promoters.bed | fold -w 60 > genes_2kb_promoters.bed.fa
以下は「macでインフィマティクス」(https://kazumaxneo.hatenablog.com/entry/2017/07/24/221538)からの引用。
flankコマンドで隣接領域を抽出している。この時-l 1000 -r 0とすることで隣接する上流1000-bp、下流0-bpが抽出されるように工夫している。ちなみにクロモソームの末端でクロモソームからはみ出すfeatureがある場合、トリミングして出力される。"
2. STREAM によるモチーフ検索
モチーフ発見ツールであるSTREAM(https://meme-suite.org/meme/tools/streme)
を用いて、モチーフ検索を行う。以下の記事が詳しいので、そちらを参照してほしい。上で作成したmulti fasta（genes_2kb_promoters.bed.fa）をそのままインプットできる。
3. GOMoによるGO解析
STREAMやMEMEと同じパイプライン（The MEME Suits）に含まれるツールGOMo（https://meme-suite.org/meme/tools/gomo）を用いて、見つけたモチーフのGenome Ontology (GO) を調べてみる。
STREAMによる解析結果のページに表示されている STREAM text output を選び、ダウンロードする。以下のようなヘッダー付きのフォーマットとなっている。
********************************************************************************
STREME - Sensitive, Thorough, Rapid, Enriched Motif Elicitation
********************************************************************************
MEME version 5.5.0 (Release date: Wed Sep 7 14:18:26 2022 -0700)
For further information on how to interpret these results please access https://meme-suite.org/meme.
To get a copy of the MEME Suite software please access https://meme-suite.org.
********************************************************************************
これをGOmoに流す。Input the motifsで先ほどの結果を選択する。今回は植物なので、SpeciesはArabidopsis thalianaを選択した。
結果はHTML形式で取得できる。
q-valueが0.05以下のGO termが得られている（閾値は変更可能）。ちなみに、 MF は Molecular Function、BP は Biological Process、CC は Cellular Component を表す。
