本手順ではgapmer型のアンチセンス核酸を例に、これが相補結合する可能性のあるRNAを検索する手順を示します (あくまでデータベースの利用例を示す為であり、何らかの保証をするものではありません)。ヒトにおけるオンターゲット、マウスにおけるターゲットの有無、そしてヒトにおけるオフターゲットを順に検索していきます。
代替手順:
RNAの配列検索において、大きく二つのアプローチが考えられます。一つは、対象となるRNA配列を準備しそれを検索するもので、直感的で理解しやすいアプローチです。しかしながら、数多くのRNAアイソフォームを対象とした検索を実施しようとすると、場合によってはゲノムの何倍もあるデータセットを構築しなければなりません (一つの遺伝子領域から、異なるエクソン組み合わせをもつRNAがたくさん転写されることもあるのです)。数多くのRNAを検索しようとすると、準備する配列データは大きくなり検索時間もかかります。もう一つのアプローチは、ゲノム配列を検索した後、RNAの座標情報(遺伝子モデル、とも呼んでいます)を使って処理をすることです。この方法では、どんなにたくさんのRNAを対象としても、配列検索を行うデータはゲノム配列よりも大きくなることはありません。
上記手順の1. ~ 3. は一つ目のアプローチ、つまりRNA配列の検索です。代替手順である4., 5. については二つ目、つまりゲノム配列の検索とその後の処理を行うアプローチに対応します。2.と4.、3.と5.はほぼ同じ結果になるはずのものですので、実際の検索ではどちらかのみで良いでしょう。ただしD3Gでは現在、アイソフォームが比較的少ないRefSeqについてはRNA配列を準備し提供していますが、他のデータセットについては、原則としてRNAの位置情報のみを提供しています。従って、D3Gのすべてのデータを使った検索を実施するのであれば、代替手順が必要になるでしょう。
ただしゲノム配列を検索するアプローチでは、一点だけ注意が必要です。RNAのエクソンはゲノム配列を鋳型として転写されたものですので、ゲノム配列の検索によってすべてのエクソンを検索できますが、エクソン・ジャンクション(エクソン同士の接合部位)配列は検索できません。こちらをも検索対象にするのであれば、やはり大きなデータセットを準備するかそれに特化した検索を補完的に実施することが必要になるでしょう。
ターミナルを開いて、検索に係るパラメータを設定します。ここではデータベースとして、RefSeqに含まれるヒトのタンパク質コード遺伝子(ミトコンドリアゲノムにコードされているものを含む)を用います。SPC–5001はLNA gapmerアンチセンス型であることから、mRNA前駆体を検索対象としています。
DB=hg38_refSeqCuratedProtCoding_prespliced_d3g1906 # GGGenomeにおける、D3G human unspliced protein-coding RNAのデータベース名
QUERY=TGCTACAAAACCCA # SPC-5001の配列。mRNAに対するアンチセンスとなっていることを銘記
STRAND=- # 検索対象とするストランド。mRNAアンチセンス鎖との一致を検索するために「-」を指定。
DISTANCE=0 # 配列検索の条件として用いる編集距離。完全一致を検索するため、「0」と指定。
TARGET=PCSK9 # ターゲット遺伝子のsymbol
OUTPUT_PREFIX=SPC-5001-PCSK9 # 出力ファイル名のプレフィックス指定パラメータを元にGGGenomeを用いた検索を行い、結果をファイルに保存します。 (検索時の条件・状況によって検索にかかる時間は異なります)
wget -O - https://gggenome.dbcls.jp/${DB}/${DISTANCE}/${STRAND}/${QUERY}.txt \
> ${OUTPUT_PREFIX}.txt検索結果が含まれるファイルの内容を表示します。
less ${OUTPUT_PREFIX}.txt以下のような内容がターミナルに表示されます。
# [ GGGenome | 2019-07-01 12:34:56 ]
# database: Human pre-spliced RNA, RefSeq curated protein coding on hg38, D3G 19.06 (Jun, 2019)
# query: TGGGTTTTGTAGCA
# count: 31
# name strand start end snippet snippet_pos snippet_end query sbjct align edit match mis del ins
NM_001042681.1|RERE|chr1:8352403:8817640:-|prespliced - 453798 453811 CCTGGAGGGGAGCTGCTTTGCCCGAGGCGGTAGCTTGCCTGTGGGCTGGGAAGAGTTCTCTGTTCTCTCCTTCGAGGGGCCCCAGGGTGGCATCCTGCTCTGGGTTTTGTAGCAAGCTCTTAGCTTCTGGCTCAGTGCAGGCTCTCCAAGGCCCCTCTTGCTCTCCAAGGCTGCTGCTTGTGGTGTCACTTCCATTGTATCAAATCTCTGCTAT 453698 453911 TGGGTTTTGTAGCA TGGGTTTTGTAGCA |||||||||||||| ============== 14 0 0 0
NM_001042682.1|RERE|chr1:8352403:8423687:-|prespliced - 59845 59858 CCTGGAGGGGAGCTGCTTTGCCCGAGGCGGTAGCTTGCCTGTGGGCTGGGAAGAGTTCTCTGTTCTCTCCTTCGAGGGGCCCCAGGGTGGCATCCTGCTCTGGGTTTTGTAGCAAGCTCTTAGCTTCTGGCTCAGTGCAGGCTCTCCAAGGCCCCTCTTGCTCTCCAAGGCTGCTGCTTGTGGTGTCACTTCCATTGTATCAAATCTCTGCTAT 59745 59958 TGGGTTTTGTAGCA TGGGTTTTGTAGCA |||||||||||||| ============== 14 0 0 0
NM_001128426.2|AP1AR|chr4:112231738:112273110:+|prespliced - 14045 14058 TACATATTTTAAAAATAAGCTTCGTTCCTGACCCATTGTAGGTCCTCAGTAAATGCAAGCCTCTCTTTGCTCTCCATTATATTAATAGCTTGTTATGTTTTGGGTTTTGTAGCATATGTAGTTGGAAACTTCAAGTCAATAAAATATTGACCTTTTTTTCAAATTGAAGAAAATAAAATGTTCTTATTATTCTCTAAACTAACAGTAATTGAGA 13945 14158 TGGGTTTTGTAGCA TGGGTTTTGTAGCA |||||||||||||| ============== 14 0 0 0
NM_001161521.1|KLHL2|chr4:165210018:165323156:+|prespliced - 85160 85173 TCTTTTGTCTCTCCACGTTTTGAATTACCTTTGTCAGCCACCTGCACTTTTCGTATAGTTGGAGAGAGAAGCTCCTTTTAGTGGGCTGTTTTTCCTTGGGTGGGTTTTGTAGCACACGGGAATCAGAGCACAAAAGCTGTGCTTAGAGGAAATGAAACTTGAGTTTCTCCTTTGTGTGCTTAGCATTTCAAGTTGTGTAAATTCTGATTCTCTT 85060 85273 TGGGTTTTGTAGCA TGGGTTTTGTAGCA |||||||||||||| ============== 14 0 0 0
NM_001161522.1|KLHL2|chr4:165210018:165323156:+|prespliced - 85160 85173 TCTTTTGTCTCTCCACGTTTTGAATTACCTTTGTCAGCCACCTGCACTTTTCGTATAGTTGGAGAGAGAAGCTCCTTTTAGTGGGCTGTTTTTCCTTGGGTGGGTTTTGTAGCACACGGGAATCAGAGCACAAAAGCTGTGCTTAGAGGAAATGAAACTTGAGTTTCTCCTTTGTGTGCTTAGCATTTCAAGTTGTGTAAATTCTGATTCTCTT 85060 85273 TGGGTTTTGTAGCA TGGGTTTTGTAGCA |||||||||||||| ============== 14 0 0 0
...以下省略...複数のエントリがヒット（検索条件に合致）していることがわかります。
これを一行ずつ確認しても構いませんが、ここでは標的遺伝子が含まれるかどうかを知りたいだけですので、
次のステップでは標的遺伝子の名前が含まれる行を抽出します。
取得した検索結果より、標的遺伝子名が含まれる行を抜き出し、その結果をファイルに保存します。
grep ${TARGET} ${OUTPUT_PREFIX}.txt \
> ${OUTPUT_PREFIX}_target.txt結果が含まれるファイルの内容を表示します。
less ${OUTPUT_PREFIX}_target.txt以下のような内容がターミナルに表示されます。
NM_174936.3|PCSK9|chr1:55039475:55064853:+|prespliced - 25311 25324 GCAAGCAGACATTTATCTTTTGGGTCTGTCCTCTCTGTTGCCTTTTTACAGCCAACTTTTCTAGACCTGTTTTGCTTTTGTAACTTGAAGATATTTATTCTGGGTTTTGTAGCATTTTTATTAATATGGTGACTTTTTAAAATAAAAACAAACAAACGTTGTCCTAAC 25211 25378 TGGGTTTTGTAGCA TGGGTTTTGTAGCA |||||||||||||| ============== 14 0 0 0結果の中に、標的遺伝子に対応するエントリが含まれていることが確認できました。
これによって、検索配列は標的遺伝子のPCSK9のmRNA前駆体の一部と完全に一致することが確認できました。
また結果が1件であることから、検索配列と一致するようなPCSK9のバリアントは検索対象であるデータベース(RefSeq)に登録されていないことがわかります。
ターミナルを開いて、検索に係る各種パラメータを設定します。ここではデータベースとして、RefSeqに含まれるマウスのタンパク質コード遺伝子(ミトコンドリアゲノムにコードされているものを含む)の前駆体を検索対象としています。
DB=mm10_refSeqCuratedProtCoding_prespliced_d3g1906 # GGGenomeにおける、D3G mouse unspliced protein-coding RNAのデータベース名。
QUERY=TGCTACAAAACCCA # SPC-5001の配列。mRNAに対するアンチセンスとなっていることを銘記。
STRAND=- # 検索対象とするストランド。mRNAアンチセンス鎖との一致を検索するために「-」を指定。
DISTANCE=2 # 配列検索の条件として用いる編集距離。ヒト・マウス間で差がある可能性を考慮し、「2」と指定。
TARGET=PCSK9 # 確認対象とする標的遺伝子名 (ヒトのオーソログであっても、マウスでは異なる名称で呼ばれている可能性もあることに注意が必要。そのようなケースかどうかを確認した上で、必要に応じてマウスの遺伝子名に修正)。
OUTPUT_PREFIX=SPC-5001-mouse # 出力ファイル名のプレフィックス。指定パラメータを元にGGGenomeを用いた検索を行い、結果をファイルに保存します。（検索時の条件・状況によって検索に要する時間は異なります）
wget -O - https://gggenome.dbcls.jp/${DB}/${DISTANCE}/${STRAND}/${QUERY}.txt \
> ${OUTPUT_PREFIX}.txt検索結果が含まれるファイルの内容を表示します。
less ${OUTPUT_PREFIX}.txtその結果、以下のような内容がターミナルに表示されます。
# [ GGGenome | 2019-07-01 12:34:56 ]
# database: Mouse pre-spliced RNA, RefSeq curated protein coding on mm10, D3G 19.06 (Jun, 2019)
# query: TGGGTTTTGTAGCA
# count: 32090
# name strand start end snippet snippet_pos snippet_end query sbjct align edit match mis del ins
NM_001001144.3|Scap|chr9:110333355:110384949:+|prespliced - 616 628 GCAACATGTAGATTCTCTGCTGGGAGTGCTCGGTTTCTGGTTTTATGTGTTGCTTAAACCTAAACTTTGCTCGCTTATGTTGCTTATGTTTACTTCCAGGTGGTTTTGTTGCATCAGGGTTATAATGGGATAGCGCTAGCTGCAGCGAGTCCCGGGTCCACACTGCATCTCAGAAAGTTGTGCAGATGGTGGCGATGAAGAGGAGGGCAGTGA 516 728 TGGGTTTTGTAGCA T-GGTTTTGTTGCA | |||||||| ||| =I========X=== 12 1 0 1
NM_001001144.3|Scap|chr9:110333355:110384949:+|prespliced - 21882 21896 CAGTGGACTGCTGGTGGTGTCTAATTTTACAGTGACACACTGCCAAGCCTGTAGACCTGTGGCTCCTTCACCTGTGTAGTACAGGATGCCCCCTAGGGAGTGGCTTTTGTAGTCACTGTTTCCTGGGTGCTTTACCATGTAGGGATTTGGAGGTAGTCTACTTTGTCTGAGCCATAATGTAGAAATGCCCTTTCTTCTTTCCTTTCTCTCTGACA 21782 21996 TGGGTTTTGTAG-CA TGGCTTTTGTAGTCA ||| |||||||| || ===X========D== 13 1 1 0
NM_001001144.3|Scap|chr9:110333355:110384949:+|prespliced - 41989 42003 TGGCCTCGAACTCAGAAATCCTCCTGCCTCTGCCTCCCAAGTGCTGGGATTAAAGGCGTGCGCCACCACCGCCTGGCCTATAAGCTTTTGTTTGTTTGTTTGGGTTTTTTTAGCAGCAAGCACTTTTAACTACTGAACTGTCTCCCCAGCCTGCACACTGCTCGCCCAGGGCTTGTGTCCAGGGAGTGCCCGATAAAGAAGCATATGTGTCTTCC 41889 42103 TGGGTTTTG-TAGCA TGGGTTTTTTTAGCA |||||||| ||||| ========XD===== 13 1 1 0
NM_001001160.3|Fbxo41|chr6:85469573:85502994:-|prespliced - 20921 20934 GTTTGAGCAGGTTGGTGGGTCTGAGCCACCATGGTTCTTGGTGTATTTCAGAGCATGAGTGTGCCTGAGAGTGAATATGAGTGCATACATAGGAACGTGTTGGGCTTTGTAGGATATGTGTTCATTGCTATGAGGCTTTCCATGTGGCTGGACATGTTCCTTTGTGAGGAGACAGTGTATTGAGGCCACAGGGAACTGGAAAAGTCTGCTTTTT 20821 21034 TGGGTTTTGTAGCA TGGGCTTTGTAGGA |||| ||||||| | ====X=======X= 12 2 0 0
NM_001001177.2|BC051142|chr17:34398819:34460734:+|prespliced - 41401 41415 TAAGGCAATTCATCACAGAATGAAAAGAAGCCTCCAGTCTTACATGTTCCCAGTGTTGTGCAATTGGAGGGCTGCCTATGACTAGCAATACAAGCCACCCTGGGTTTTCCTAGCATTTCCATGCACAATGTTGAAAGATGGGTATTCAGTTATTAGGTATTACCTTCACATCAATTAAGCATAATGATCTCTCTCCCCATCTCTCTAGCCGGAGC 41301 41515 TGGGTTTTG-TAGCA TGGGTTTTCCTAGCA |||||||| ||||| ========XD===== 13 1 1 0
...以下省略...取得した検索結果より、標的遺伝子名が含まれる行を抜き出し、その結果をファイルに保存します。
grep -i \|${TARGET}\| ${OUTPUT_PREFIX}.txt \
> ${OUTPUT_PREFIX}_target.txt抜き出された結果を表示します。
less ${OUTPUT_PREFIX}_target.txtその結果、以下のような内容がターミナルに表示されます。
<出力なし>結果が得られていないことから、mouseのPCSK9前駆体mRNA配列には、2塩基以内の違いで一致する領域が無いことがわかりました。
ターミナルを開いて、検索に係る各種パラメータを設定します。ここでは手順1–1.と同様に、データベースとしてRefSeqに含まれるヒトのタンパク質コード遺伝子(ミトコンドリアゲノムにコードされているものを含む)を用います。
DB=hg38_refSeqCuratedProtCoding_prespliced_d3g1906 # GGGenomeにおける、D3G human unspliced protein-coding RNAのデータベース名
QUERY=TGCTACAAAACCCA # SPC-5001の配列。mRNAに対するアンチセンスとなっていることを銘記
STRAND=- # 検索対象とするストランド。mRNAアンチセンス鎖との一致を検索するために「-」を指定。
DISTANCE=2 # 配列検索の条件として用いる編集距離。オフターゲットの影響も考えるため、ここでは「2」と指定。
OUTPUT_PREFIX=SPC-5001-prespliced # 出力ファイル名のプレフィックス指定パラメータを元にGGGenomeを用いた検索を行い、結果をファイルに保存します。（検索時の条件・状況によって検索に要する時間は異なります）
wget -O - https://gggenome.dbcls.jp/${DB}/${DISTANCE}/${STRAND}/${QUERY}.txt \
> ${OUTPUT_PREFIX}.txt検索結果が含まれるファイルの内容を表示します。
less ${OUTPUT_PREFIX}.txt次のような内容がターミナルに表示されます。
# [ GGGenome | 2019-07-01 12:34:56 ]
# database: Human pre-spliced RNA, RefSeq curated protein coding on hg38, D3G 19.06 (Jun, 2019)
# query: TGGGTTTTGTAGCA
# count: 66949
# name strand start end snippet snippet_pos snippet_end query sbjct align edit match mis del ins
NM_000019.3|ACAT1|chr11:108121530:108148168:+|prespliced - 19646 19658 AAATTAGCTGGGTATAGCAGCACACGCCTGTAGTCCCAGCTACTAGGGAATCCAAGGCAGGAGGATCAATTGAGCCCAGGAGTCCGAGGCTGCAGTGAGCTAGGTTTGTAGCATTGCACTCTATTGGGAGACTGTGCAAGACCATTTCAAACAAACAAAATTAGCTGGCTGTGGTGGCTCCCACCTGTAGTCCCAGCTACTTGGGAGGCTGAG 19546 19758 TGGGTTTTGTAGCA TAGG-TTTGTAGCA | || ||||||||| =X==I========= 12 1 0 1
NM_000020.2|ACVRL1|chr12:51907417:51923361:+|prespliced - 15624 15637 CTAGGTTCCCAGATCATTAGGGCAGAGTTTGCACGTCCTCTGGTCACTGGAATCCACCCAGCCCACGAATCATCTCCCTCTTGAAGGATTTTATTTCTACTGGGTTTTGGAACAAACTCCTGCTGAGACCCCACAGCCAGAAACTGAAAGCAGCAGCTCCCCAAAGCCTGGAAAATCCCTAAGAGAAGGCCTGGGGCAGGAAGTGGAGTGACAG 15524 15737 TGGGTTTTGTAGCA TGGGTTTTGGAACA ||||||||| | || =========X=X== 12 2 0 0
NM_000022.3|ADA|chr20:44619518:44651758:-|prespliced - 5170 5182 GTACCTCAGGTAGACAACTCCTGAAACTCAGCTTCCCCCTCTGTAAAATGGGGTGACAAAACCAAGATCTTGGGGTTCTTGGGGAAACTGACATGCTGATTGGTTTTTGTACAGTGCCTGGCTGGTAACAGCAGGCCCTCAGGGGTGCGTTTCCTTCCTGGGGACTGGAGTGGGGGTTGCAGTAGACTCTGGGAGGCCTCTCCAGCTGCAGAA 5070 5282 TGGGTTTTGTAGCA TGGTTTTTGTA-CA ||| ||||||| || ===X=======I== 12 1 0 1
NM_000033.3|ABCD1|chrX:153724867:153744762:+|prespliced - 7174 7187 CACCGTCTCTATCTGATGGTTTGGTTAGCTTCAGTTCTTGGGAGTATCATTCGCCTGTCTTCTGTCTGCTGGTCTTCACTCATGGTGGATCATTTCCTTGTGGGTTTTGCAACAGAACTTCTATGGGGATTGGTGTGGCCTCAATGGAGGGTGTAAGTCCAAATGTCCTTCTGTTTTGGCCAAGCACTCCAGGGTAAGCAGTTGGAGGCCATTA 7074 7287 TGGGTTTTGTAGCA TGGGTTTTGCAACA ||||||||| | || =========X=X== 12 2 0 0
NM_000034.3|ALDOA|chr16:30053089:30070420:+|prespliced - 7583 7595 TTATTAATGACAAGATCTAAACCCTTCCAAGGGAGAGAGCCAAGAAATGCTTTTAATACCCTCTTTTCTGAGGTAATGGCAGCAGAGCCTCTCTCAGCCATGGGTTCTGTGCACAGAATGATGGTGTGGTCCTTCAAAAGCCCGCAGAAGTGATTTCATCCCCAGACTGAGACTCCCTAGAGACACCCCAAATCACACCCTCATCAAACACAA 7483 7695 TGGGTTTTGTAGCA TGGGTTCTGT-GCA |||||| ||| ||| ======X===I=== 12 1 0 1
...以下省略...一致配列が多数存在することがわかります。ただ、同一遺伝子であっても異なるバリアントに存在する一致配列、あるいは同一バリアントであっても異なる領域に存在しうる一致配列もすべて異なる行として現れていることに注意が必要です。
取得した検索結果を遺伝子名で並び替えつつ、遺伝子毎に区切ってファイルに保存します。
tail -n +6 ${OUTPUT_PREFIX}.txt \
| cut -f 1-4,6- \
| sort -k 2,2 -t '|' \
| awk 'BEGIN{FS="|"}{if (prev!=$2){print "\n>>>"$2} prev=$2;print}' \
> ${OUTPUT_PREFIX}_per_gene.txt遺伝子単位にまとめた結果ファイルの内容を表示します。
less ${OUTPUT_PREFIX}_per_gene.txt次のような内容がターミナルに表示されます。
>>>A1CF
NM_001198818.1|A1CF|chr10:50799408:50885675:-|prespliced - 10821 10833 10721 10933 TGGGTTTTGTAGCA TGGTTTTTGT-GCA ||| |||||| ||| ===X======I=== 12 1 0 1
NM_001198818.1|A1CF|chr10:50799408:50885675:-|prespliced - 63557 63569 63457 63669 TGGGTTTTGTAGCA TGGGTTTTCTA-CA |||||||| || || ========X==I== 12 1 0 1
NM_001198818.1|A1CF|chr10:50799408:50885675:-|prespliced - 83779 83792 83679 83892 TGGGTTTTGTAGCA TTGGTTTTGTAGGA | |||||||||| | =X==========X= 12 2 0 0
NM_001198819.1|A1CF|chr10:50799408:50885675:-|prespliced - 10821 10833 10721 10933 TGGGTTTTGTAGCA TGGTTTTTGT-GCA ||| |||||| ||| ===X======I=== 12 1 0 1
NM_001198819.1|A1CF|chr10:50799408:50885675:-|prespliced - 63557 63569 63457 63669 TGGGTTTTGTAGCA TGGGTTTTCTA-CA |||||||| || || ========X==I== 12 1 0 1
NM_001198819.1|A1CF|chr10:50799408:50885675:-|prespliced - 83779 83792 83679 83892 TGGGTTTTGTAGCA TTGGTTTTGTAGGA | |||||||||| | =X==========X= 12 2 0 0
NM_001198820.1|A1CF|chr10:50799408:50885675:-|prespliced - 10821 10833 10721 10933 TGGGTTTTGTAGCA TGGTTTTTGT-GCA ||| |||||| ||| ===X======I=== 12 1 0 1
NM_001198820.1|A1CF|chr10:50799408:50885675:-|prespliced - 63557 63569 63457 63669 TGGGTTTTGTAGCA TGGGTTTTCTA-CA |||||||| || || ========X==I== 12 1 0 1
NM_001198820.1|A1CF|chr10:50799408:50885675:-|prespliced - 83779 83792 83679 83892 TGGGTTTTGTAGCA TTGGTTTTGTAGGA | |||||||||| | =X==========X= 12 2 0 0
NM_014576.3|A1CF|chr10:50799408:50885675:-|prespliced - 10821 10833 10721 10933 TGGGTTTTGTAGCA TGGTTTTTGT-GCA ||| |||||| ||| ===X======I=== 12 1 0 1
NM_014576.3|A1CF|chr10:50799408:50885675:-|prespliced - 63557 63569 63457 63669 TGGGTTTTGTAGCA TGGGTTTTCTA-CA |||||||| || || ========X==I== 12 1 0 1
NM_014576.3|A1CF|chr10:50799408:50885675:-|prespliced - 83779 83792 83679 83892 TGGGTTTTGTAGCA TTGGTTTTGTAGGA | |||||||||| | =X==========X= 12 2 0 0
NM_138932.2|A1CF|chr10:50799408:50885675:-|prespliced - 10821 10833 10721 10933 TGGGTTTTGTAGCA TGGTTTTTGT-GCA ||| |||||| ||| ===X======I=== 12 1 0 1
NM_138932.2|A1CF|chr10:50799408:50885675:-|prespliced - 63557 63569 63457 63669 TGGGTTTTGTAGCA TGGGTTTTCTA-CA |||||||| || || ========X==I== 12 1 0 1
NM_138932.2|A1CF|chr10:50799408:50885675:-|prespliced - 83779 83792 83679 83892 TGGGTTTTGTAGCA TTGGTTTTGTAGGA | |||||||||| | =X==========X= 12 2 0 0
NM_138933.2|A1CF|chr10:50799408:50885675:-|prespliced - 10821 10833 10721 10933 TGGGTTTTGTAGCA TGGTTTTTGT-GCA ||| |||||| ||| ===X======I=== 12 1 0 1
NM_138933.2|A1CF|chr10:50799408:50885675:-|prespliced - 63557 63569 63457 63669 TGGGTTTTGTAGCA TGGGTTTTCTA-CA |||||||| || || ========X==I== 12 1 0 1
NM_138933.2|A1CF|chr10:50799408:50885675:-|prespliced - 83779 83792 83679 83892 TGGGTTTTGTAGCA TTGGTTTTGTAGGA | |||||||||| | =X==========X= 12 2 0 0
>>>A2ML1
NM_144670.5|A2ML1|chr12:8822471:8876783:+|prespliced - 10167 10180 10067 10280 TGGGTTTTGTAGCA TGGGTTTTATAACA |||||||| || || ========X==X== 12 2 0 0
>>>A4GNT
NM_016161.2|A4GNT|chr3:138123717:138132387:-|prespliced - 4644 4656 4544 4756 TGGGTTTTGTAGCA TGGGTTTTCT-GCA |||||||| | ||| ========X=I=== 12 1 0 1
...以下省略...検索結果が遺伝子単位にまとまっていることがわかります。
参考までに、いくつの遺伝子にヒットしているのか調べます。まずは、遺伝子名だけの行を抽出します。
grep '>>>' ${OUTPUT_PREFIX}_per_gene.txt | sort > ${OUTPUT_PREFIX}_per_gene.genes.txtその後、行を数えます。
cat ${OUTPUT_PREFIX}_per_gene.genes.txt | wc -l結果が次のように表示されます。
75767,000個以上の遺伝子にヒットしていることがわかりました。
ターミナルを開いて、検索に係る各種パラメータを設定します。検索対象として、リファレンス・ヒトゲノムを用います。検索により得られた配列一致が、遺伝子領域に存在するかどうかを判断するため、遺伝子のゲノム座標データについても準備します。
wget http://d3g.riken.jp/release/latest/hg38/ncbiRefSeqCuratedProteinCoding.bed12.gz # 遺伝子のゲノム座標データの準備
GENE_COORDINATES=ncbiRefSeqCuratedProteinCoding.bed12.gz # 遺伝子のゲノム座標データファイル
DB=hg38 # GGGenomeにおける、ヒトhg38ゲノムのデータベース名
QUERY=TGCTACAAAACCCA # SPC-5001の配列。mRNAに対するアンチセンスとなっていることを銘記
STRAND= # 検索対象とするストランド。ゲノム上には両鎖に遺伝子が存在し得ることから、片鎖に限定せず両鎖を検索する。
DISTANCE=0 # 配列検索の条件として用いる編集距離。完全一致を検索するため、「0」と指定。
TARGET=PCSK9 # ターゲット遺伝子のsymbol
OUTPUT_PREFIX=SPC-5001-PCSK9 # 出力ファイル名のプレフィックス指定パラメータを元にGGGenomeを用いた検索を行い、結果をファイルに保存します。 (検索時の条件・状況によって検索にかかる時間は異なります)
wget -O - https://gggenome.dbcls.jp/${DB}/${DISTANCE}/${STRAND}/${QUERY}.bed \
> ${OUTPUT_PREFIX}.bed検索結果が含まれるファイルの内容を表示します。
less ${OUTPUT_PREFIX}.bed以下のような内容がターミナルに表示されます。
track name=GGGenome description="GGGenome matches"
chr1 8363829 8363843 . 0 +
chr1 177713585 177713599 . 0 +
chr1 184489057 184489071 . 0 +
chr3 27217664 27217678 . 0 +
chr3 55813946 55813960 . 0 +
chr4 150736744 150736758 . 0 +
chr5 44364503 44364517 . 0 +
chr7 153934740 153934754 . 0 +
chr12 31992431 31992445 . 0 +
...以下省略...ゲノム上に、多数の配列一致が見つかったことがわかります。
配列一致の見つかったゲノム位置が、遺伝子と重複する領域に存在するかを調べ、結果をファイルに保存します。
bedtools intersect -a ${OUTPUT_PREFIX}.bed -b ${GENE_COORDINATES} -wa -wb -S \
> ${OUTPUT_PREFIX}_gene.bed結果が含まれるファイルの内容を表示します。
less ${OUTPUT_PREFIX}_gene.bed次のような内容がターミナルに表示されます。
chr1 8363829 8363843 . 0 + chr1 8352403 8423687 NM_001042682.1|RERE 0 - 8355086 8364133 0,0,0 13 2717,181,147,721,223,1379,114,162,200,93,163,81,127, 0,3015,3696,5792,7360,7708,9359,10279,11652,12342,13408,70323,71157,
chr1 8363829 8363843 . 0 + chr1 8352403 8817640 NM_012102.3|RERE 0 - 8355086 8656297 0,0,0 24 2717,181,147,721,223,1379,114,162,200,93,163,81,99,100,125,49,105,97,106,126,71,469,185,481, 0,3015,3696,5792,7360,7708,9359,10279,11652,12342,13408,70323,113521,142659,145001,156223,188810,204071,205014,262157,271906,303569,440000,464756,
chr1 8363829 8363843 . 0 + chr1 8352403 8817640 NM_001042681.1|RERE 0 - 8355086 8656297 0,0,0 23 2717,181,147,721,223,1379,114,162,200,93,163,81,99,100,125,49,105,97,106,126,71,469,481, 0,3015,3696,5792,7360,7708,9359,10279,11652,12342,13408,70323,113521,142659,145001,156223,188810,204071,205014,262157,271906,303569,464756,
chr3 27217664 27217678 . 0 + chr3 27110902 27369460 NM_152534.4|NEK10 0 - 27111271 27352882 0,0,0 39 418,56,53,109,111,101,30,141,38,55,87,184,214,71,130,76,103,122,46,138,129,103,65,78,62,140,225,87,80,68,79,42,85,99,131,61,108,40,236, 0,5037,5172,8857,20977,30579,32542,51558,51798,60916,63536,63747,81126,90607,91525,145393,173699,173937,176795,179714,180359,180581,182685, 184710,186276,190793,193844,196956,198023,200046,201196,203394,211274,233369,235183,241562,241909,256620,258322,
chr3 27217664 27217678 . 0 + chr3 27215605 27369460 NM_199347.3|NEK10 0 - 27215831 27352882 0,0,0 25 275,76,103,122,46,138,129,103,65,78,62,140,225,87,80,68,79,42,85,99,131,61,108,40,236, 0,40690,68996,69234,72092,75011,75656,75878,77982,80007,81573,86090,89141,92253,93320,95343,96493,98691,106571,128666,130480,136859,137206,151917,153619,
chr3 55813946 55813960 . 0 + chr3 55508310 56468363 NM_015576.2|ERC2 0 - 55683832 56435007 0,0,0 18 2966,66,135,148,161,136,12,194,141,141,138,168,168,156,75,417,797,116, 0,175483,191067,226460,380078,442114,477666,483746,498870,502138,510583,572506,631198,640666,665135,787708,926040,959937,
chr4 150736744 150736758 . 0 + chr4 150264658 151015267 NM_001199282.2|LRBA 0 - 150265721 151014642 0,0,0 57 1154,152,197,102,165,156,63,178,90,168,153,120,141,113,116,103,118,137,147,125,167,109,65,62,134,79,134,442,160,108,122,181,154,179,1059,193,124,82,109,93,98,63,80,169,153,109,134,198,147,120,127,122,96,101,232,435,66, 0,13194,17791,21274,37969,45570,50902,56532,61150,8 5333,150779,170930,172065,203014,206965,223073,226259,323389,326054,334348,418892,470599,497124,533422,541612,543661,552465,563521,567158,579441,579999,584159,584763,586065,587226,603012,603523,605866,606686,608004,628393,631735,633080,635390,641179,641638,643675,644001,649536,650949,651742,651958,656539,663857,664175,749768,750543,
chr4 150736744 150736758 . 0 + chr4 150264658 151015497 NM_006726.4|LRBA 0 - 150265721 151014642 0,0,0 58 1154,152,197,102,168,156,63,178,90,168,153,120,141,113,116,103,118,137,147,33,125,167,109,65,62,134,79,134,442,160,108,122,181,154,179,1059,193,124,82,109,93,98,63,80,169,153,109,134,198,147,120,127,122,96,101,232,435,25, 0,13194,17791,21274,37966,45570,50902,56532,61150,8 5333,150779,170930,172065,203014,206965,223073,226259,323389,326054,332424,334348,418892,470599,497124,533422,541612,543661,552465,563521,567158,579441,579999,584159,584763,586065,587226,603012,603523,605866,606686,608004,628393,631735,633080,635390,641179,641638,643675,644001,649536,650949,651742,651958,656539,663857,664175,749768,750814,
chr5 44364503 44364517 . 0 + chr5 44304994 44388682 NM_004465.1|FGF10 0 - 44304994 44388682 0,0,0 3 198,104,325, 0,5432,83363,
chr1 55064785 55064799 . 0 - chr1 55039475 55064853 NM_174936.3|PCSK9 0 + 55039837 55063584 0,0,0 12 569,192,124,134,142,197,184,174,149,178,182,1485, 0,4367,7047,12802,13174,16517,17855,18560,19023,20010,21899,23893,
...以下省略...重複を調べた結果から、標的配列遺伝子名を含む行を抜き出し、その結果をファイルに保存します。
grep ${TARGET} ${OUTPUT_PREFIX}_gene.bed \
> ${OUTPUT_PREFIX}_gene_target.txt抜き出された結果を表示します。
less ${OUTPUT_PREFIX}_gene_target.txtその結果、以下のような内容がターミナルに表示されます。
chr1 55064785 55064799 . 0 - chr1 55039475 55064853 NM_174936.3|PCSK9 0 + 55039837 55063584 0,0,0 12 569,192,124,134,142,197,184,174,149,178,182,1485, 0,4367,7047,12802,13174,16517,17855,18560,19023,20010,21899,23893,ヒトゲノム上で検索配列と一致した領域群に、標的遺伝子のPCSK9がコードされている領域が含まれていることが確認できました。
ターミナルを開いて、検索に係る各種パラメータを設定します。手順4–1.と同様に、検索対象としてリファレンス・ヒトゲノムを用います。検索により得られた配列一致が、遺伝子領域に存在するかどうかを判断するため、遺伝子のゲノム座標データについても準備します(手順4–1.で準備したものと同じものですので、そちらの実施後でしたら下のwgetコマンドは必要ありません)。
wget -nc http://d3g.riken.jp/release/latest/hg38/ncbiRefSeqCuratedProteinCoding.bed12.gz # 遺伝子座標データの準備
GENE_COORDINATES=ncbiRefSeqCuratedProteinCoding.bed12.gz # 遺伝子のゲノム座標データファイル
DB=hg38 # GGGenomeにおける、ヒトhg38ゲノムのデータベース名
QUERY=TGCTACAAAACCCA # SPC-5001の配列。mRNAに対するアンチセンスとなっていることを銘記
STRAND= # 検索対象とするストランド。ゲノム上には両鎖に遺伝子が存在し得ることから、片鎖に限定せず両鎖を検索する。
DISTANCE=2 # 配列検索の条件として用いる編集距離。配列検索の条件として用いる編集距離。オフターゲットの影響も考えるため、ここでは「2」と指定。
OUTPUT_PREFIX=SPC-5001-prespliced # 出力ファイル名のプレフィックスwget -O - https://gggenome.dbcls.jp/${DB}/${DISTANCE}/${STRAND}/${QUERY}.bed \
> ${OUTPUT_PREFIX}.bed結果ファイルの内容を表示します。
less ${OUTPUT_PREFIX}.bedその結果、以下のような内容がターミナルに表示されます。
track name=GGGenome description="GGGenome matches"
chr1 12428 12442 . 0 +
chr1 21471 21484 . 0 +
chr1 56174 56189 . 0 +
chr1 173165 173177 . 0 +
chr1 182947 182961 . 0 +
chr1 191995 192008 . 0 +
chr1 460993 461006 . 0 +
chr1 674286 674299 . 0 +
chr1 695964 695977 . 0 +
chr1 807983 807995 . 0 +
... 以下省略 ...bedtoolsを用いてゲノム上の検索結果と、遺伝子領域とが重複する領域だけを抜き出し、その結果をファイルに保存します。
bedtools intersect -a ${OUTPUT_PREFIX}.bed -b ${GENE_COORDINATES} -wa -wb -S -f 0.85 \
> ${OUTPUT_PREFIX}_gene_intersect.txtアンチセンス配列をクエリとしているため-Sオプションを指定しています。また、デフォルトの場合だと一塩基でも重複している場合も含めて出力されますので、一致配列の85%以上が重複するもの（つまり、14塩基長のうち12塩基以上が重複）のみを出力するよう「-f 0.85」とオプションを指定しています。
また、重複塩基長がオフターゲット候補にならないほど短いと考えられるケース（ゲノム上の遺伝子領域の末端から外にはみ出してマッチしているケース）を除くため、マッチ領域の50%以上が遺伝子領域と重複する箇所のみが抽出されるように-f 0.5のオプションを指定しています。
結果が含まれるファイルの内容を表示します。
less ${OUTPUT_PREFIX}_gene_intersect.txt次のような内容がターミナルに表示されます。
chr1 1072070 1072082 . 0 + chr1 1071745 1074307 NM_001205252.1|RNF223 0 - 1071816 1072566 0,0,0 2 830,292, 0,2270,
chr1 1391261 1391274 . 0 + chr1 1385710 1398342 NM_001350497.1|CCNL2 0 - 1387230 1390858 0,0,0 10 1872,93,112,142,105,100,125,65,144,110, 0,2066,2243,4519,4748,5055,6968,7685,9683,12522,
chr1 1391261 1391274 . 0 + chr1 1385710 1399342 NM_001350499.1|CCNL2 0 - 1387230 1391562 0,0,0 13 1872,93,112,142,105,100,279,125,65,121,110,75,324, 0,2066,2243,4519,4748,5055,5586,6968,7685,9683,12522,12886,13308,
chr1 1391261 1391274 . 0 + chr1 1385710 1399342 NM_001350498.1|CCNL2 0 - 1387230 1390858 0,0,0 12 1872,93,112,142,105,100,125,65,121,110,75,324, 0,2066,2243,4519,4748,5055,6968,7685,9683,12522,12886,13308,
chr1 1391261 1391274 . 0 + chr1 1385710 1399342 NM_001320153.2|CCNL2 0 - 1387230 1390858 0,0,0 12 1872,93,112,142,105,810,14,65,121,110,75,324, 0,2066,2243,4519,4748,5055,7079,7685,9683,12522,12886,13308,
chr1 1391261 1391274 . 0 + chr1 1385710 1399342 NM_001350500.1|CCNL2 0 - 1387230 1391562 0,0,0 13 1872,93,112,142,105,100,279,14,65,121,110,75,324, 0,2066,2243,4519,4748,5055,5586,7079,7685,9683,12522,12886,13308,
chr1 1391261 1391274 . 0 + chr1 1385710 1399342 NM_001320155.2|CCNL2 0 - 1387230 1390858 0,0,0 12 1872,93,112,142,105,100,14,65,121,110,75,324, 0,2066,2243,4519,4748,5055,7079,7685,9683,12522,12886,13308,
chr1 1391261 1391274 . 0 + chr1 1385710 1399342 NM_030937.5|CCNL2 0 - 1387230 1399306 0,0,0 11 1872,93,112,142,105,100,65,121,110,75,324, 0,2066,2243,4519,4748,5055,7685,9683,12522,12886,13308,
chr1 1545828 1545842 . 0 + chr1 1541672 1574882 NM_014188.2|SSU72 0 - 1542065 1574557 0,0,0 5 495,119,140,144,405, 0,2196,3190,23100,32805,
chr1 1648125 1648139 . 0 + chr1 1635225 1659097 NM_033489.2|CDK11B 0 - 1635763 1657264 0,0,0 21 631,190,149,117,108,122,106,122,91,176,66,159,96,123,137,139,128,116,47,124,184, 0,711,1107,1456,1671,1855,2182,2536,3274,5051,5822,6436,6680,7160,9900,14273,17213,20143,22030,22149,23688,
... 以下省略 ...重複する遺伝子名を行頭に出力した後、これを辞書順に並べます（コマンドの１〜５行目）。
これによって同じ遺伝子に対してゲノム上で重複する一致配列が連続するようになります。
その後、異なる遺伝子が現れる毎に区切りを出力し（コマンドの６行目）、その結果をファイルに保存します。
cat ${OUTPUT_PREFIX}_gene_intersect.txt \
| awk '{print $10"\t"$0}' \
| cut -f 2- -d '|' \
| sort -k1 -k2 -k3 \
| awk '{if (prev!=$1){print "\n>>>"$1} prev=$1;print}' \
> ${OUTPUT_PREFIX}_per_gene_intersect.txt結果ファイルの内容を表示します。
less ${OUTPUT_PREFIX}_per_gene_intersect.txtその結果、以下のような内容がターミナルに表示されます。
>>>A1CF
A1CF chr10 50801883 50801897 . 0 + chr10 50799408 50885675 NM_001198818.1|A1CF 0 - 50806728 50859940 0,0,0 14 7472,149,137,182,274,98,165,239,131,135,144,93,48,95, 0,10485,11631,14448,16597,21143,28722,36665,42453,44579,60433,63480,64624,86172,
A1CF chr10 50801883 50801897 . 0 + chr10 50799408 50885675 NM_001198819.1|A1CF 0 - 50806728 50850787 0,0,0 15 7472,149,137,206,274,98,165,239,131,135,143,144,93,48,95, 0,10485,11631,14448,16597,21143,28722,36665,42453,44579,51256,60433,63480,64624,86172,
A1CF chr10 50801883 50801897 . 0 + chr10 50799408 50885675 NM_001198820.1|A1CF 0 - 50806728 50850787 0,0,0 14 7472,149,137,182,274,98,165,239,131,135,143,144,48,95, 0,10485,11631,14448,16597,21143,28722,36665,42453,44579,51256,60433,64624,86172,
A1CF chr10 50801883 50801897 . 0 + chr10 50799408 50885675 NM_014576.3|A1CF 0 - 50806728 50859940 0,0,0 13 7472,149,137,182,274,98,165,239,131,135,144,48,95, 0,10485,11631,14448,16597,21143,28722,36665,42453,44579,60433,64624,86172,
A1CF chr10 50801883 50801897 . 0 + chr10 50799408 50885675 NM_138932.2|A1CF 0 - 50806728 50859940 0,0,0 13 7472,149,137,206,274,98,165,239,131,135,144,48,95, 0,10485,11631,14448,16597,21143,28722,36665,42453,44579,60433,64624,86172,
A1CF chr10 50801883 50801897 . 0 + chr10 50799408 50885675 NM_138933.2|A1CF 0 - 50806728 50850787 0,0,0 13 7472,149,137,182,274,98,165,239,131,135,143,144,95, 0,10485,11631,14448,16597,21143,28722,36665,42453,44579,51256,60433,86172,
A1CF chr10 50822106 50822119 . 0 + chr10 50799408 50885675 NM_001198818.1|A1CF 0 - 50806728 50859940 0,0,0 14 7472,149,137,182,274,98,165,239,131,135,144,93,48,95, 0,10485,11631,14448,16597,21143,28722,36665,42453,44579,60433,63480,64624,86172,
A1CF chr10 50822106 50822119 . 0 + chr10 50799408 50885675 NM_001198819.1|A1CF 0 - 50806728 50850787 0,0,0 15 7472,149,137,206,274,98,165,239,131,135,143,144,93,48,95, 0,10485,11631,14448,16597,21143,28722,36665,42453,44579,51256,60433,63480,64624,86172,
A1CF chr10 50822106 50822119 . 0 + chr10 50799408 50885675 NM_001198820.1|A1CF 0 - 50806728 50850787 0,0,0 14 7472,149,137,182,274,98,165,239,131,135,143,144,48,95, 0,10485,11631,14448,16597,21143,28722,36665,42453,44579,51256,60433,64624,86172,
A1CF chr10 50822106 50822119 . 0 + chr10 50799408 50885675 NM_014576.3|A1CF 0 - 50806728 50859940 0,0,0 13 7472,149,137,182,274,98,165,239,131,135,144,48,95, 0,10485,11631,14448,16597,21143,28722,36665,42453,44579,60433,64624,86172,
A1CF chr10 50822106 50822119 . 0 + chr10 50799408 50885675 NM_138932.2|A1CF 0 - 50806728 50859940 0,0,0 13 7472,149,137,206,274,98,165,239,131,135,144,48,95, 0,10485,11631,14448,16597,21143,28722,36665,42453,44579,60433,64624,86172,
A1CF chr10 50822106 50822119 . 0 + chr10 50799408 50885675 NM_138933.2|A1CF 0 - 50806728 50850787 0,0,0 13 7472,149,137,182,274,98,165,239,131,135,143,144,95, 0,10485,11631,14448,16597,21143,28722,36665,42453,44579,51256,60433,86172,
A1CF chr10 50874842 50874855 . 0 + chr10 50799408 50885675 NM_001198818.1|A1CF 0 - 50806728 50859940 0,0,0 14 7472,149,137,182,274,98,165,239,131,135,144,93,48,95, 0,10485,11631,14448,16597,21143,28722,36665,42453,44579,60433,63480,64624,86172,
A1CF chr10 50874842 50874855 . 0 + chr10 50799408 50885675 NM_001198819.1|A1CF 0 - 50806728 50850787 0,0,0 15 7472,149,137,206,274,98,165,239,131,135,143,144,93,48,95, 0,10485,11631,14448,16597,21143,28722,36665,42453,44579,51256,60433,63480,64624,86172,
A1CF chr10 50874842 50874855 . 0 + chr10 50799408 50885675 NM_001198820.1|A1CF 0 - 50806728 50850787 0,0,0 14 7472,149,137,182,274,98,165,239,131,135,143,144,48,95, 0,10485,11631,14448,16597,21143,28722,36665,42453,44579,51256,60433,64624,86172,
A1CF chr10 50874842 50874855 . 0 + chr10 50799408 50885675 NM_014576.3|A1CF 0 - 50806728 50859940 0,0,0 13 7472,149,137,182,274,98,165,239,131,135,144,48,95, 0,10485,11631,14448,16597,21143,28722,36665,42453,44579,60433,64624,86172,
A1CF chr10 50874842 50874855 . 0 + chr10 50799408 50885675 NM_138932.2|A1CF 0 - 50806728 50859940 0,0,0 13 7472,149,137,206,274,98,165,239,131,135,144,48,95, 0,10485,11631,14448,16597,21143,28722,36665,42453,44579,60433,64624,86172,
A1CF chr10 50874842 50874855 . 0 + chr10 50799408 50885675 NM_138933.2|A1CF 0 - 50806728 50850787 0,0,0 13 7472,149,137,182,274,98,165,239,131,135,143,144,95, 0,10485,11631,14448,16597,21143,28722,36665,42453,44579,51256,60433,86172,
>>>A2ML1
A2ML1 chr12 8832637 8832651 . 0 - chr12 8822471 8876783 NM_144670.5|A2ML1 0 + 8822651 8875011 0,0,0 36 242,184,163,53,21,160,85,127,115,110,168,228,61,146,150,195,91,115,229,127,122,52,84,177,82,157,75,163,215,216,128,91,69,103,42,726, 0,710,1248,7255,12190,13035,13783,14968,15864,16641,18897,20662,22970,23605,25077,26248,27197,27688,29312,29738,31656,32308,33037,34692,35035,35474,38409,38663,41322,45370,45758,46065,46663,51953,52499,53586,
>>>A4GNT
A4GNT chr3 138127731 138127744 . 0 + chr3 138123717 138132387 NM_016161.2|A4GNT 0 - 138124263 138131256 0,0,0 3 1161,434,176, 0,7131,8494,
... 以下省略...ゲノムでなく、タンパク質コードRNA前駆体の配列を、直接検索したときと同じ配列が見つかっていることが確認できます (１１カラム目)。また、２-４カラム目はクエリ配列がヒットしたゲノム座標を示していますが、A1CFに対する上から6つのエントリはいずれも同じ値を示しています。これは、ゲノム上の同じ箇所に一致したものであること、そしてその領域にはRNA前駆体の構造が複数登録されていることを示しています。
3–3と同様に、いくつの遺伝子にヒットしているのか調べます。まずは、遺伝子名だけの行を抽出します。
grep '>>>' ${OUTPUT_PREFIX}_per_gene_intersect.txt| sort > ${OUTPUT_PREFIX}_per_gene_intersect.genes.txtその後、行を数えます。
cat ${OUTPUT_PREFIX}_per_gene_intersect.genes.txt | wc -l結果が次のように表示されます。
75757,000個以上の遺伝子にヒットしていることがわかりました。
ただ、3–5では7576個の遺伝子であったにもかかわらず、一つ足りないようです。どこが違うのかを確認してみます。
diff ${OUTPUT_PREFIX}_per_gene.genes.txt ${OUTPUT_PREFIX}_per_gene_intersect.genes.txtすると次のように結果が表示されます。
624d623
< >>>BAGE5BAGE5、という遺伝子が手順3–5 (ゲノムでなくmRNA前駆体を対象に検索したもの) のみに見つかっていることがわかりました。
これがどんなヒットであったのか、確認してみます。
grep "BAGE5" ${OUTPUT_PREFIX}_per_gene.txt次のような結果が出力されます。
>>>BAGE5
NM_182484.2|BAGE5|chr13_KN538372v1_fix:76209:170143:+|prespliced - 21817 21830 AATGAGATACATTTATATTTTTTTGCATGGTAACTTATAAGGCATTGCCAAGAGAGGGAATTCATTGGAGGCAAAGAAAAATATAAAACATAGTGTATATAGGGCTTTGTAGCATCCGTGGTTTCAGGCATTTATGGGAGCCAGGCAGGGGTGGGGTATGGGGTTGAAATGTATTTCCTATGGATAAGCACAGGACTACTTTATACCACTTAAT 21717 21930 TGGGTTTTGTAGCA AGGGCTTTGTAGCA ||| ||||||||| X===X========= 12 0chr13_KN538372v1_fix、という見慣れない染色体の上に定義されていることがわかりました。これを個別に調べてみると、どうやらリファレンス・ゲノム配列の小規模な訂正（パッチ、hg38 patch2 ）として、リファレンス・ゲノムの発表後に追加された配列であることがわかりました。一方、GGGenomeを用いて実施したhg38に対する検索は、パッチによる訂正を含まない配列に対するものでした。以上から、二つのアプローチによる検索は基本的に等価であること、ただ遺伝子セットが定義されるタイミングなどによって結果は僅かに異なるケースがあることがわかりました。