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

いったんファイルをチェックしてみる。

Sbicolor_454_v3.1.1.fa.fai

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

中身を確認。

genome.txt

Chr01	80884392
Chr02	77742459
Chr03	74386277
Chr04	68658214
Chr05	71854669
Chr06	61277060
Chr07	65505356
Chr08	62686529
Chr09	59416394
Chr10	61233695

次に、ターゲット遺伝子のリストを作成する。スタートはgff3ファイルとした。Phytozomeなどから入手できます。

Sbicolor_454_v3.1.1.genes.gff

##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) などを想定しています。

my_genes.txt

Sobic.005G047100
Sobic.003G270300
Sobic.001G192800
Sobic.010G101900
Sobic.002G374700
Sobic.005G035100

次に、以下のスクリプトを用いて情報を抽出・整理します。

make_bed_from_gff3.py

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は、こんな感じ。

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=に一致させる必要がある。

my_genes_2.txt

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を作っている。

my_genes_2.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 を表す。