[[ノート/ノート]]~
訪問者数 &counter();      最終更新 &lastmod();

[[ノート/breseq/チュートリアルを試す1]] > [[ノート/breseq/チュートリアルを試す2]] > [[ノート/breseq/UT-Austinの進化データを試す]]~
[[ノート/breseq/Multi-Isolate1]] > [[ノート/breseq/Multi-Isolate2]]

*breseq Multi-Isolate その1 2017-04-01 [#udbbd631]

DDBJから、multi-isolateのプロジェクトを探して、breseq処理をしてみる

出発点は、multi-isolateだと変異の同定が難しいのではないか、と言う疑問。

** DDBJのBioProjectのプロパティでmonoisolateでないものを探す [#iead4c17]

DDBJのデータベースの関係は、[[BioProject, BioSample, DRA へのデータ登録:https://www.genome-sci.jp/old2010-2015/pdf/20140820furuya.pdf]]の7枚目のスライドに図示されている。

この中の、BioProjectデータベースには、プロジェクトを登録されているが、
そのプロパティの中に、Project typeとして、Sample Scope が Monoisolate であるという情報がある。欲しいのはヘテロな環境でのNGSデータなので、Monoisolateではないものを探す。

データベースの記述によると、Multiisolateというのがあるらしい。この種類のデータを実験材料にしたい。

それで、プロジェクトを1つずつチェックしていくのも手ではあるが、面倒なので、
データベースのデータをダウンロードして、ローカルで検索することにする。

データは、
-[[BioProject(プロジェクト一覧)のホーム:http://trace.ddbj.nig.ac.jp/bioproject/index.html]]から、上部タブの「Download」をクリック
-[[/ddbj_database/bioproject のインデックス:ftp://ftp.ddbj.nig.ac.jp/ddbj_database/bioproject]] に飛ぶが、その中でまずは[[ddbj_core_bioproject.xml:ftp://ftp.ddbj.nig.ac.jp/ddbj_database/bioproject/ddbj_core_bioproject.xml]]をダウンロードして検索する。もし全体で探したければbioproject.xmlがよいのだろう(大きい)。

XMLなのだが、面倒ならテキストファイルとしてエディタで検索してもよい。
       <ProjectType>
         <ProjectTypeSubmission>
           <Target capture="eWhole" material="eGenome" sample_scope="eMonoisolate">
             <Organism species="139" taxID="224326">
               <OrganismName>Borrelia burgdorferi B31</OrganismName>
というような顔をしているので、
ねらい目は<Target>の中のsample_scopeがeMultiisolateであるかどうか、である。
また、更に解析対象物を知りたいので、同じ<Target>の中の<OrganismName>を取り出しておく。

pythonにxmlパーザーがあるので、ちょっとしたプログラムを書けば済む。たとえば、
 #!/usr/bin/env python
 # coding: UTF-8
 import xml.etree.ElementTree as ET
 #tree = ET.parse("bioproject.xml")
 tree = ET.parse("ddbj_core_bioproject.xml")
 root = tree.getroot()
 for u in root.iter('Package'):
     str = ''
     for v in u.iter('ArchiveID'):
         str += v.get('accession')
     for v in u.iter('Target'):
         str += ',  ' + v.get('sample_scope')
         #text = ''
         for w in v.iter('OrganismName'):
             text = (w.text).strip()
         str += ',  ' + text
     print(str )

この出力をファイルに取ったうえで、sample_scopeのところがMultiisolateである行
だけをgrepなどで拾った。プログラム内でフィルタしてもいいのだが、全体が見たかった
ので2段階処理にしてみた。

結果は、&ref(acc_multi.list,,このよう);になった。


ここから、プロジェクトを選ぼうと思う。
登録されているプロジェクトの1つに、PRJDB3087があって、
eMultiisolateであり、 Escherichia coliを対象としているので、おもしろそうだ。

------------------------------------------------------
**PRJDB3087にあるデータをマップしてみる [#h3e0b1d2]

[[DRA Search:http://trace.ddbj.nig.ac.jp/DRASearch/]]により、キーワードに PRJDB3087 を入れて検索すると[[3つの結果:http://trace.ddbj.nig.ac.jp/DRASearch/query?keyword=PRJDB3087&show=20]]が得られ、その中の
DRA002503.experiment.xml を見ると、3つのAccession 
DRX021052 、DRX021053 、DRX021054 がある。
それぞれのリンクをクリックしてページを開くと、
-[[DRX021052:http://trace.ddbj.nig.ac.jp/DRASearch/experiment?acc=DRX021052]]からFASTQとして[[DRR022997_1.fastq.bz2:ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/DRA002/DRA002503/DRX021052/DRR022997_1.fastq.bz2]]、[[DRR022997_2.fastq.bz2:ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/DRA002/DRA002503/DRX021052/DRR022997_2.fastq.bz2]] が
-[[DRX021053:http://trace.ddbj.nig.ac.jp/DRASearch/experiment?acc=DRX021053]]からFASTQとして[[DRR022998_1.fastq.bz2:ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/DRA002/DRA002503/DRX021053/DRR022998_1.fastq.bz2]]、[[DRR022998_2.fastq.bz2:ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/DRA002/DRA002503/DRX021053/DRR022998_2.fastq.bz2]] が
-[[DRX021054:http://trace.ddbj.nig.ac.jp/DRASearch/experiment?acc=DRX021054]]からFASTQとして[[DRR022999_1.fastq.bz2:ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/DRA002/DRA002503/DRX021054/DRR022999_1.fastq.bz2]]、[[DRR022999_2.fastq.bz2:ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/DRA002/DRA002503/DRX021054/DRR022999_2.fastq.bz2]] が

が得られる。

これらの6つのデータそれぞれについてクオリティチェックをした結果と、2つずつのペアリードデータを併せてbreseqでマップした結果を、以下に置く。手順は
 fastqc *.bz2
 bunzip2 DRR022997_1.fastq.b2z 他のファイルも同じ
 breseq -j 16 -o output -r REL606.gbk DRR022997_1.fastq DRR022997_2.fastq
breseqによるマッピング処理について、予め以下の点に疑問がある。
-クオリティによるトリミングについて、検討していないので、マッピングがうまく行かない、もしくは品質が悪くなる可能性がある(高い)。
-リファレンスゲノムはよくわからないのでbreseqの元論文と同じREL606を使ってしまったが、当然「外れ」が多くなるであろうから、もっと調べる方がよかろう。とにかく試してみたレベルである。

|DRR022997 | [[DRR022997_1 FastQC:http://ginger.yy.is.sci.toho-u.ac.jp/~yamanouc/NGS2017/DRR022997_1_fastqc.html]] | [[DRR022997_2 FastQC:http://ginger.yy.is.sci.toho-u.ac.jp/~yamanouc/NGS2017/DRR022997_2_fastqc.html]] |[[breseqの結果:http://ginger.yy.is.sci.toho-u.ac.jp/~yamanouc/NGS2017/DRR022997_output]] |
|DRR022998 | [[DRR022998_1 FastQC:http://ginger.yy.is.sci.toho-u.ac.jp/~yamanouc/NGS2017/DRR022998_1_fastqc.html]] | [[DRR022998_2 FastQC:http://ginger.yy.is.sci.toho-u.ac.jp/~yamanouc/NGS2017/DRR022998_2_fastqc.html]] | [[breseqの結果:http://ginger.yy.is.sci.toho-u.ac.jp/~yamanouc/NGS2017/DRR022998_output]] |
|DRR022999 | [[DRR022999_1 FastQC:http://ginger.yy.is.sci.toho-u.ac.jp/~yamanouc/NGS2017/DRR022999_1_fastqc.html]] | [[DRR022999_2 FastQC:http://ginger.yy.is.sci.toho-u.ac.jp/~yamanouc/NGS2017/DRR022999_2_fastqc.html]] | [[breseqの結果:http://ginger.yy.is.sci.toho-u.ac.jp/~yamanouc/NGS2017/DRR022999_output]] |

***Trimmingをかけてみる [#u8346cc5]
Trimmingの効果を見るため、trimmaticで下記の処理をする
-クオリティ20未満のリードを除去
-長さ36未満のリードを除去
-ウィンドウサイズ4で見た時その中の平均クオリティが15未満のリードを除去

 DRR22997の結果
 TrimmomaticPE: Started with arguments:
  -threads 16 -phred33 ../DRR022997_1.fastq ../DRR022997_2.fastq
   DRR022997_1_trim20.fastq DRR022997_1_unpaired.fastq
   DRR022997_2_trim20.fastq DRR022997_2_unpaired.fastq
   TRAILING:20 MINLEN:36 SLIDINGWINDOW:4:15
 Input Read Pairs: 2496778 Both Surviving: 2456608 (98.39%)
   Forward Only Surviving: 34113 (1.37%) Reverse Only Surviving: 4714 (0.19%)
   Dropped: 1343 (0.05%)
 TrimmomaticPE: Completed successfully

 DRR22998の結果
 TrimmomaticPE: Started with arguments:
  -threads 16 -phred33 ../DRR022998_1.fastq ../DRR022998_2.fastq
   DRR022998_1_trim20.fastq DRR022998_1_unpaired.fastq
   DRR022998_2_trim20.fastq DRR022998_2_unpaired.fastq
   TRAILING:20 MINLEN:36 SLIDINGWINDOW:4:15
 Input Read Pairs: 1151142 Both Surviving: 1137953 (98.85%)
   Forward Only Surviving: 12861 (1.12%) Reverse Only Surviving: 137 (0.01%)
   Dropped: 191 (0.02%)
 TrimmomaticPE: Completed successfully

 DRR22999の結果
 TrimmomaticPE: Started with arguments:
  -threads 16 -phred33 ../DRR022999_1.fastq ../DRR022999_2.fastq
   DRR022999_1_trim20.fastq DRR022999_1_unpaired.fastq
   DRR022999_2_trim20.fastq DRR022999_2_unpaired.fastq
   TRAILING:20 MINLEN:36 SLIDINGWINDOW:4:15
 Input Read Pairs: 1262796 Both Surviving: 1249828 (98.97%)
   Forward Only Surviving: 12617 (1.00%) Reverse Only Surviving: 171 (0.01%)
   Dropped: 180 (0.01%)
 TrimmomaticPE: Completed successfully

改めてFastQCでクオリティをチェック。上記ログのDroppedがそれほど多くない割には、FastQCの形はまあきれいになっていると言えるか。頭(左)の側の乱れはまだかなり大きいと感じる。~
このトリムしたデータを、breseqでマップしてみる。

|DRR022997_trim20 | [[DRR022997_1 FastQC:http://ginger.yy.is.sci.toho-u.ac.jp/~yamanouc/NGS2017/DRR022997_1_trim20_fastqc.html]] | [[DRR022997_2 FastQC:http://ginger.yy.is.sci.toho-u.ac.jp/~yamanouc/NGS2017/DRR022997_2_trim20_fastqc.html]] |[[breseqの結果:http://ginger.yy.is.sci.toho-u.ac.jp/~yamanouc/NGS2017/DRR022997_output]] |
|DRR022998_trim20 | [[DRR022998_1 FastQC:http://ginger.yy.is.sci.toho-u.ac.jp/~yamanouc/NGS2017/DRR022998_1_trim20_fastqc.html]] | [[DRR022998_2 FastQC:http://ginger.yy.is.sci.toho-u.ac.jp/~yamanouc/NGS2017/DRR022998_2_trim20_fastqc.html]] | [[breseqの結果:http://ginger.yy.is.sci.toho-u.ac.jp/~yamanouc/NGS2017/DRR022998_output]] |
|DRR022999_trim20 | [[DRR022999_1 FastQC:http://ginger.yy.is.sci.toho-u.ac.jp/~yamanouc/NGS2017/DRR022999_1_trim20_fastqc.html]] | [[DRR022999_2 FastQC:http://ginger.yy.is.sci.toho-u.ac.jp/~yamanouc/NGS2017/DRR022999_2_trim20_fastqc.html]] | [[breseqの結果:http://ginger.yy.is.sci.toho-u.ac.jp/~yamanouc/NGS2017/DRR022999_output]] |


------------------------

***メモ [#q0af1980]

[[ncbiのBioProject内でST131をサーチ:https://www.ncbi.nlm.nih.gov/bioproject/?term=ST131]]  36エントリ。かなりたくさんmultiisolateがある。

以下、いくつか目を引いたもの。

[[Global diversity and evolutionary history of Escherichia coli sequence type 131 (ST131):https://www.ncbi.nlm.nih.gov/bioproject/297860]]  これも面白いかな?~
表のNumber of Linksのエントリー部分をクリックすると、データへも行けるみたい。~
[[Links from BioProject:https://www.ncbi.nlm.nih.gov/sra?linkname=bioproject_sra_all&from_uid=297860]] 138このプロジェクトデータ~
[[SRX1460182: Whole genome sequencing of E. coli ST131:https://www.ncbi.nlm.nih.gov/sra/SRX1460182(accn) ]]  1つ目のデータ例~
[[SRX1460182: Whole genome sequencing of E. coli ST131:https://www.ncbi.nlm.nih.gov/sra/SRX1460182]]  1つ目のデータ例~
[[SRR2970783:https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR2970783]]~
SRA-Toolkit(のfastq-dump)を使えとさ。


[[PRJDB4303:https://www.ncbi.nlm.nih.gov/bioproject/336158]]  Comparative genomics of ESBL-producing Escherichia coli ST131 isolates

[[PRJDB3868:https://www.ncbi.nlm.nih.gov/bioproject/330822]]  Whole genome sequencing of Escherichia coli ST131

[[PRJNA311313:https://www.ncbi.nlm.nih.gov/bioproject/311313]]  Escherichia coli JJ1887 Genome sequencing
[[PRJNA307507:https://www.ncbi.nlm.nih.gov/bioproject/307507]]  Escherichia coli Genome sequencing and assembly Comparative sequencing of ExPEC ST131

トップ   編集 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS