ノート/breseq/Multi-Isolate１

サイトトップ

ノート
データ処理
├マイニング
├バスケット解析
├バスケット解析をＲで
├図書貸出をＲで
├モデルクラスタリング
└
テキストマイニング
├テキストマイニングTM
├TMとシソーラス
├PubMedをTM
├TMとMeSH
├TMとNLTK
├テキストマイニングとtagger
├医薬品添付文書ＤＢ
├論文の処理1
├剽窃１
├Bing検索
└
分子進化学
├分子進化学
└
Pythonと論文アクセス
├Pythonを使ってみる
├PythonでPubMed
├Pythonで...続き
└
Macメモ
├Script

&ref(): File not found: "favicon_story.txt" at page "ノート/ノート";

自信ないので
├Winのバッチ
│├スクリプト入門
│├バッチファイル
│├SETコマンド
│├自分のメモ
│└
└

メニュー編集

breseq Multi-Isolate その１　2017-04-01 †

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

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

↑

DDBJのBioProjectのプロパティでmonoisolateでないものを探す †

DDBJのデータベースの関係は、BioProject, BioSample, DRA へのデータ登録の７枚目のスライドに図示されている。

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

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

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

データは、

BioProject（プロジェクト一覧）のホームから、上部タブの「Download」をクリック
/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などで拾った。プログラム内でフィルタしてもいいのだが、全体が見たかったので２段階処理にしてみた。

結果は、このようになった。

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

↑

PRJDB3087にあるデータをマップしてみる †

DRA Searchにより、キーワードに PRJDB3087 を入れて検索すると３つの結果が得られ、その中の DRA002503.experiment.xml を見ると、３つのAccession　 DRX021052 、DRX021053 、DRX021054　がある。それぞれのリンクをクリックしてページを開くと、

DRX021052からFASTQとしてDRR022997_1.fastq.bz2、DRR022997_2.fastq.bz2　が
DRX021053からFASTQとしてDRR022998_1.fastq.bz2、DRR022998_2.fastq.bz2　が
DRX021054からFASTQとしてDRR022999_1.fastq.bz2、DRR022999_2.fastq.bz2　が

が得られる。

これらの６つのデータそれぞれについてクオリティチェックをした結果と、２つずつのペアリードデータを併せて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	DRR022997_2 FastQC	breseqの結果
DRR022998	DRR022998_1 FastQC	DRR022998_2 FastQC	breseqの結果
DRR022999	DRR022999_1 FastQC	DRR022999_2 FastQC	breseqの結果

↑

Trimmingをかけてみる †

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	DRR022997_2 FastQC	breseqの結果
DRR022998_trim20	DRR022998_1 FastQC	DRR022998_2 FastQC	breseqの結果
DRR022999_trim20	DRR022999_1 FastQC	DRR022999_2 FastQC	breseqの結果

↑

メモ †

ncbiのBioProject内でST131をサーチ　36エントリ。かなりたくさんmultiisolateがある。

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

Global diversity and evolutionary history of Escherichia coli sequence type 131 (ST131)　　これも面白いかな？
表のNumber of Linksのエントリー部分をクリックすると、データへも行けるみたい。
Links from BioProject　138このプロジェクトデータ
SRX1460182: Whole genome sequencing of E. coli ST131　　１つ目のデータ例
SRR2970783
SRA-Toolkit（のfastq-dump）を使えとさ。

PRJDB4303　　Comparative genomics of ESBL-producing Escherichia coli ST131 isolates

PRJDB3868　　Whole genome sequencing of Escherichia coli ST131

PRJNA311313　　Escherichia coli JJ1887 Genome sequencing PRJNA307507　　Escherichia coli Genome sequencing and assembly　Comparative sequencing of ExPEC ST131

最新の5件

breseq Multi-Isolate その１ 2017-04-01 †

DDBJのBioProjectのプロパティでmonoisolateでないものを探す †

PRJDB3087にあるデータをマップしてみる †

Trimmingをかけてみる †

メモ †

breseq Multi-Isolate その１　2017-04-01 †