Pythonバイオ?　Pythonバイオ/ツール?
6047　　　2019-05-12 (日) 16:31:26

リードデータ取得 †

ここでは、既に公開されている実験データを使ってマッピング処理を試すことを考えている。手元で実験を行いシーケンサーからの出力データがある場合は、当然ながらそのデータを使うことになる。

A comprehensive comparison of RNA-Seq-based transcriptome analysis from reads to differential gene expression and cross-comparison with microarrays: a case study in Saccharomyces cerevisiae　Saccharomycesの例を使った例

GEO Accession viewer GSE37599

↑

リードデータを検索したい場合 †

↑

例�　〔臈帖屮肇薀鵐好�リプトーム解析」　p71 でのリードデータ取得 †

この例の場合は、NBCIの公開実験データベースSRA（Sequence Read Archive）から投稿全体（SRA、submission accession）のアクセション番号であるSRA番号 SRAをアクセスする（p72）

興味のあるデータが、submission_accession=SRA000299であるとする。（論文等から与えられている）

これから、NCBIのSRAデータベースでSRA000299を検索する。

Webで手作業でアクセスする場合は、直接にNCBIのSRAデータベースにアクセスすればよい。

https://www.ncbi.nlm.nih.gov/sra/?term=SRA000299

以下、Web-APIベース（＝Pythonプログラム経由）で作業することを考える。NCBIの SRAデータベースはWeb-API（のやり方）を公開していないようで、Entrezを使うように 指示している。（Download SRA sequences from Entrez search results）　 ここの例にあるように、Entrezの問い合わせを作る。

Entrez APIのメモ

Download SRA sequences from Entrez search results

SRA Toolkit download <--fastq-dump and sam-dump

PythonでNCBIのAPIから文献情報を取得してみた ではBioPythonを使わないで直接Web-APIを叩いているが、今はこれは使わないでおく。

PythonによるEntrez SRAのアクセス

from Bio import Entrez
Entrez.email = "yamanouc@hyperresearch.com"
handle = Entrez.esearch(db="sra", term="SRA000299")
result = Entrez.read(handle)    # 検索結果をresultに入れる
print('SRA SRA000299, IDList', result['IdList'])
for id in result['IdList'][:1]:    # resultのうちIdListを取り出して、１つずつアクセス
    handle = Entrez.efetch(db="sra", id=id, retmode="xml")  # １つずつアクセス
    print(handle.read())

結果はxmlしかない。xmlを解読して、欲しいRUN accession id を入手するには、

import xml.etree.ElementTree as ET
from Bio import Entrez

Entrez.email = "yamanouc@hyperresearch.com"
handle = Entrez.esearch(db="sra", term="SRA000299")
result = Entrez.read(handle)
print('SRA SRA000299, IDList', result['IdList'])
for id in result['IdList']:
    handle = Entrez.efetch(db="sra", id=id, retmode="xml")
    root = ET.fromstring(handle.read())
    # find
    items = root.findall('.//RUN')
    for i, u in enumerate(items):
        #print(i, 'tag', u.tag, 'attr', u.attrib)
        print(u.get('accession'))

結果は

SRR002324
SRR002320
SRR002325
SRR002322
SRR002321
SRR002323

となる。

このrun accession number SRR002320〜5を使って、fastqテーブルをアクセスする。

ここから、SRAのデータ（SRAデータ、最後にはfastqデータ）のダウンロード。
まずSRAのページ⇒NCBI SRA Toolkit
ここからNCBI SRA Toolkitをダウンロード・展開して、これによってアクセスする。Rだと中で自動的にいろいろとやってくれるようだ。

fasterq-dump SRR002320

結果は

spots read      : 39,266,713
reads read      : 39,266,713
reads written   : 39,266,713

で、ファイルは SRR002320.fastqという8340468900バイトのファイルが作られた。先頭を覗いてみると

@SRR002320.1 080226_CMLIVERKIDNEY_0007:1:1:112:735 length=36
GTGGTGGGGTTGGTATTTGGTTTCTCGTTTTAATTA
+SRR002320.1 080226_CMLIVERKIDNEY_0007:1:1:112:735 length=36
IIIIIIII"IIIII)I$I1%HII"I#./(#/'$#*#
@SRR002320.2 080226_CMLIVERKIDNEY_0007:1:1:114:564 length=36
GGATACTCAGGCTGGCCCAATTTCTGGGCGTGGGAA
+SRR002320.2 080226_CMLIVERKIDNEY_0007:1:1:114:564 length=36
IIII:>&<I;I%I88II1&+I:IF>II,&D:I-'),
@SRR002320.3 080226_CMLIVERKIDNEY_0007:1:1:109:558 length=36
GTAGAATTAGAATTGTGAAGATGATAAGTGTAGAGG
+SRR002320.3 080226_CMLIVERKIDNEY_0007:1:1:109:558 length=36
IIIIIIIIIIIIIIIIIIIIIII<IIAIIII6I?I:

のようになっている。

リファレンス配列の方は、NCBI GenomeにおいてHomo sapiens genomeを検索
　⇒ 検索
　　⇒ GRCh38.p12 (December 2017) Download　⇒ ファイルGRCh38.p12.tarとして作る(938MB)　　(p90) 　　　

↑

遺伝研講習会での例 †

「先進ゲノム支援」情報解析講習会のご案内
⇒情報解析講習会ビデオ＜2018年度　情報解析講習会（中級者向け）＞
⇒資料（GitHub）

1-1の題材は出芽酵母Saccharomyces cerevisiaeで2つの異なる条件での培養
Intawat Nookaew et al
"A comprehensive comparison of RNA-Seq-based transcriptome analysis from reads to differential gene expression and cross-comparison with microarrays: a case study in Sacchaomyces cerevisiae"
Nucleic Acids Research, 2012, Vol. 40, No. 20, Septemter 2012 doi: 10.1093/nar/gks804
full text

論文での記述（p10095）に、リードに関するACCESSION NUMBERSとして

GSE37599, SRS307298, SRR453566, SRR453567,
SRR453568, SRR453569, SRR453570, SRR453571 and
SRR453578.

とあるので、これを頼りに、上の例と同じようにして、SRAのデータ（SRAデータ、最後にはfastqデータ）をダウンロードできる。 まずSRAのページ ⇒ NCBI SRA Toolkit からNCBI SRA Toolkitをダウンロード・展開して、これを使ってアクセスする。SRR453566のデータのダウンロードは

fasterq-dump SRR453566
のようにすればよい。

また論文中に referenceとしてS288cゲノムを使うことが書かれているので、このデータも アクセスする必要がある。ここを見た。

Transcriptome analysis using reference genome-based reads mapping The genome sequence of S. cerevisiae strain S288c and its annotations were retrieved from the SGD databases and used for all analysis.

ここから、SGD (Saccharomyces Genome Database) のS288Cをサーチすると、Strain: S288Cが得られる。この中で、GenBank GCF_000146045.2のエントリーを使うことにする。

GCF_000146045.2をクリックすると、GenBankのR64 Organism name: Saccharomyces cerevisiae S288C (baker's yeast) Strain: S288Cのページに跳ぶ。跳んだ先のページの右側「Access the data」のDownload the RefSeq assemblyをクリックすると、ファイルリストが表示される。その中から、参照シーケンスとしてGCF_000146045.2_R64_genomic.fna.gz、記述GFFファイルとしてGCF_000146045.2_R64_genomic.gff.gzをダウンロードする。

別の到達パスとしては、NBCIのgenomeのページからsaccharomyces cerevisiae s288c[orgn] を検索する。
この中のFastaフォーマットのgenomeをダウンロードをクリックすると、ファイルGCF_000146045.2_R64_genomic.fna.gzが得られる。
更にgffフォーマットのannotationをダウンロードをクリックすると、ファイルGCF_000146045.2_R64_genomic.gff.gzが得られる。

↑

その他いろいろ †

RNA-seq演習(2018-03)高橋弘喜

↑

Ensembl Genomesのアクセス †

↑

UCSC ? †

↑

FastQCによるクォリティチェック †

バッチモードで実行する場合。

fastqc --nogroup SRR453566_1.fastq
fastqc --nogroup SRR453566_2.fastq

クオリティのグラフが生成される。使い方の詳細はFASTQ クオリティコントロールにパラメータ指定あり。

元ページはfastqc

↑

Trimmomaticによるトリミング †

Trimmomaticはこちら参照。
ダウンロード（2019-03-02時点でVersion 0.38）
/usr/local/Trimmomaticにインストール

使い方は

java -jar trimmomatic-0.38.jar PE -phred33 input_forward.fq.gz input_reverse.fq.gz
output_forward_paired.fq.gz output_forward_unpaired.fq.gz 
output_reverse_paired.fq.gz output_reverse_unpaired.fq.gz 
ILLUMINACLIP:TruSeq3- PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

This will perform the following:

Remove adapters (ILLUMINACLIP:TruSeq3-PE.fa:2:30:10)
Remove leading low quality or N bases (below quality 3) (LEADING:3)
Remove trailing low quality or N bases (below quality 3) (TRAILING:3)
Scan the read with a 4-base wide sliding window, cutting when the average quality
  per base drops below 15 (SLIDINGWINDOW:4:15)
Drop reads below the 36 bases long (MINLEN:36)

遺伝研スライドによると

java -jar -Xmx512m trimmomatic-0.38.jar \ 
  PE                     \ 
  -threads ${NSLOTS}         \ 
  -phred33               \ 
  -trimlog log_SRR${NUM}.txt \ 
  SRR${NUM}_1.fastq.gz  \ 
  SRR${NUM}_2.fastq.gz  \ 
  paired_SRR${NUM}_1.trim.fastq.gz   \ 
  unpaired_SRR${NUM}_1.trim.fastq.gz \ 
  paired_SRR${NUM}_2.trim.fastq.gz   \ 
  unpaired_SRR${NUM}_2.trim.fastq.gz \ 
  ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10 \ 
  LEADING:20 \ 
  TRAILING:20 \ 
  SLIDINGWINDOW:4:15 \ 
  MINLEN:36

一般に、アダプタ除去は他の加工より先に行う方が良い（他の加工によりマッチングが難しくなるため）としている。

個々の細かいパラメータの意味はtrimmomaticのページとマニュアルに書かれている。以下、上記の例について説明する（Paired Endの場合に限る）。

-Xmx256mこれはjavaコマンドに対する指定で、trimmomaticのパラメータではない。-Xmx256mはメモリ割当ての最大量を256Mバイトにする。無指定時は64M
PE動作モードがSE(SingleEnd)かPE(PairedEnd)か
--threads 16処理時の並列スレッド数
-phred33塩基（読取り）品質の記述法、-phread33か-phread64、無指定時は自動判別（-v0.32以降）
-trimlog log_SRR453566.txt実行ログの出力先ファイル名の指定
SRR453566_1.fastqPairedEndでの入力forwardファイル
SRR453566_2.fastqPairedEndでの入力backwardファイル
paired_SRR453566}_1.trim.fastqPairedEndでの出力paired forwardファイル
unpaired_SRR453566}_1.trim.fastqPairedEndでの出力unpaired forwardファイル
paired_SRR453566}_2.trim.fastqPairedEndでの出力paired backwardファイル
unpaired_SRR453566}_2.trim.fastqPairedEndでの出力unpaired backwardファイル

これ以降は、個別の除去ステップを指定する。

ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10ステップ1でIllumina adapterを除去、TruSeq***はアダプターを記述したfastaファイル、2は最大ミスマッチ数、30は回文アライメント時に２つの隣接リードがどれだけ正確にマッチするかを指定、10はアダプターとリード間のアライメントマッチの正確さ
LEADING:20先頭から低品質ベースを取り除く、この時の残すための最低品質が20
TRAILING:20末尾から低品質ベースを取り除く、この時残すための最低品質が20
CROP:? 例では使われていない品質に関係なく、先頭から指定された塩基数だけを残し後ろを除去
HEADCROP:? 例では使われていない品質に関係なく、先頭から指定された塩基数だけ除去し後ろを残す
SLIDINGWINDOW:4:15スライディングウィンドウ幅を4とし、その中での平均品質が15以上のものを残す
MINLEN:36（通常最後に行う）残っているリードのうち、長さの最小値36以上のものを残す

実際のコマンドは

java -jar -Xmx512m /usr/local/Trimmomatic/trimmomatic-0.38.jar PE \
 -threads 32         \
 -phred33               \
 -trimlog log_SRR453566.txt \
 SRR453566_1.fastq  \
 SRR453566_2.fastq  \
 paired_SRR453566_1.trim.fastq   \
 unpaired_SRR453566_1.trim.fastq \
 paired_SRR453566_2.trim.fastq   \
 unpaired_SRR453566_2.trim.fastq \
 ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10 \
 LEADING:20 \
 TRAILING:20 \
 SLIDINGWINDOW:4:15 \
 MINLEN:36

結果は

TrimmomaticPE: Started with arguments:
 -threads 32 -phred33 -trimlog log_SRR453566.txt SRR453566_1.fastq SRR453566_2.fastq 
paired_SRR453566}_1.trim.fastq unpaired_SRR453566_1.trim.fastq 
paired_SRR453566_2.trim.fastq unpaired_SRR53566_2.trim.fastq ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10
LEADING:20 TRAILING:20 SLIDINGWINDOW:4:15 MINLEN:36
Using PrefixPair: 'TACACTCTTTCCCTACACGACGCTCTTCCGATCT' and 
'GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT'
Using Long Clipping Sequence: 'AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA'
Using Long Clipping Sequence: 'AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC'
Using Long Clipping Sequence: 'GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT'
Using Long Clipping Sequence: 'TACACTCTTTCCCTACACGACGCTCTTCCGATCT'
ILLUMINACLIP: Using 1 prefix pairs, 4 forward/reverse sequences, 0 forward only 
sequences, 0 reverse only sequences
Input Read Pairs: 5725730 Both Surviving: 5115482 (89.34%) Forward Only Surviving: 
514793 (8.99%) Reverse Only Surviving: 46123 (0.81%) Dropped: 49332 (0.86%)
TrimmomaticPE: Completed successfully

なお、IlluminaClipで指定するアダプターシーケンスは GitHubのtrimmomaticのパッケージ中のadapters/TruSeq30PE-2.fa を使うことができた。

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2019-05-12 (日) 16:31:26 (1417d)