[[Pythonバイオ]] [[Pythonバイオ/ツール]]~
&counter();   &lastmod();~

*リードデータ取得 [#a9d72733]

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

[[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:https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3488244/]] Saccharomycesの例を使った例

[[GEO Accession viewer GSE37599:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE37599]]



**リードデータを検索したい場合 [#i763f2bc]
-[[公共NGSデータの検索と登録(2017NGSハンズオン講習会-8月29日):https://biosciencedbc.jp/gadget/human/20170829_3_nakazato_20170818.pdf]]

-[[DDBJ Sequence Read Archive Handbook:https://www.ddbj.nig.ac.jp/dra/submission.html#metadata]]  [[データベースはどこ?:https://yokazaki.hatenablog.com/entry/2015/05/27/215757]]

-[[門田先生 農学生命情報科学 特論I 第1回 データベース、データ取得、ファイル形式、Quality Control:http://www.iu.a.u-tokyo.ac.jp/~kadota/20150616_kadota.pdf]]
**例 〔臈帖屮肇薀鵐好リプトーム解析」 p71 でのリードデータ取得 [#j1940e2a]

この例の場合は、NBCIの公開実験データベースSRA(Sequence Read Archive)から投稿全体(SRA、submission accession)のアクセション番号であるSRA番号
SRAをアクセスする(p72)~
// [[Getting data from the SRA:https://edwards.sdsu.edu/research/getting-data-from-the-sra/]]

興味のあるデータが、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):https://www.ncbi.nlm.nih.gov/sra/docs/sradownload/]] 
ここの例にあるように、Entrezの問い合わせを作る。

Entrez APIのメモ
-[[NCBI APIs:https://www.ncbi.nlm.nih.gov/home/develop/api/]]
-[[Entrez Programming Utilities Help:https://www.ncbi.nlm.nih.gov/books/NBK25501/]]
-[[The E-utilities In-Depth: Parameters, Syntax and More:https://www.ncbi.nlm.nih.gov/books/NBK25499/#chapter4.EFetch]]
-[[Table 1 – Valid values of &retmode and &rettype for EFetch (null = empty string):https://www.ncbi.nlm.nih.gov/books/NBK25499/table/chapter4.T._valid_values_of__retmode_and/?report=objectonly]]

[[Download SRA sequences from Entrez search results:https://www.ncbi.nlm.nih.gov/sra/docs/sradownload/]]

[[SRA Toolkit download:https://www.ncbi.nlm.nih.gov/sra/docs/toolkitsoft/]] <--fastq-dump and sam-dump

[[PythonでNCBIのAPIから文献情報を取得してみた:https://qiita.com/MTNakata/items/1538da3e97fe0a8b951a]] では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を取り出して、1つずつアクセス
     handle = Entrez.efetch(db="sra", id=id, retmode="xml")  # 1つずつアクセス
     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のページ:https://www.ncbi.nlm.nih.gov/sra]]⇒[[NCBI SRA Toolkit:https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=software]]~
ここから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を検索~
 ⇒ [[検索:https://www.ncbi.nlm.nih.gov/search/all/?term=Homo+sapiens+genome]]~
  ⇒ GRCh38.p12 (December 2017) Download ⇒ ファイルGRCh38.p12.tarとして作る(938MB)  (p90)
   


*遺伝研講習会での例 [#c85616f7]
[[「先進ゲノム支援」情報解析講習会のご案内:https://www.genome-sci.jp/whatsnew/event/news20180920.html]]~
⇒[[情報解析講習会ビデオ<2018年度 情報解析講習会(中級者向け)>:https://www.genome-sci.jp/lecture20181st]]~
⇒[[資料(GitHub):https://github.com/genome-sci/python_bioinfo_2018]]

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:https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3488244/]]

論文での記述(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:https://www.yeastgenome.org/]] (Saccharomyces Genome Database) のS288Cをサーチすると、[[Strain: S288C:https://www.yeastgenome.org/strain/S000203483]]が得られる。この中で、GenBank GCF_000146045.2のエントリーを使うことにする。

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


別の到達パスとしては、NBCIのgenomeの[[ページ:https://www.ncbi.nlm.nih.gov/genome]]からsaccharomyces cerevisiae s288c[orgn] を[[検索:https://www.ncbi.nlm.nih.gov/genome?term=saccharomyces+cerevisiae+s288c%5Borgn%5D&cmd=DetailsSearch]]する。~
この中の[[Fastaフォーマットのgenomeをダウンロード:ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/146/045/GCF_000146045.2_R64/GCF_000146045.2_R64_genomic.fna.gz]]をクリックすると、ファイルGCF_000146045.2_R64_genomic.fna.gzが得られる。~
更に[[gffフォーマットのannotationをダウンロード:ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/146/045/GCF_000146045.2_R64/GCF_000146045.2_R64_genomic.gff.gz]]をクリックすると、ファイルGCF_000146045.2_R64_genomic.gff.gzが得られる。




* その他いろいろ [#v944cf06]

[[RNA-seq演習:https://www.genome-sci.jp/seminar2018/7_doc_takahashi.pdf]](2018-03)高橋弘喜

**Ensembl Genomesのアクセス [#l037d90e]

**UCSC ? [#d97a4303]


*FastQCによるクォリティチェック [#k7363eae]
バッチモードで実行する場合。
 fastqc --nogroup SRR453566_1.fastq
 fastqc --nogroup SRR453566_2.fastq

クオリティのグラフが生成される。
クオリティのグラフが生成される。使い方の詳細は[[FASTQ クオリティコントロール:https://bi.biopapyrus.jp/rnaseq/qc/fastqc.html]]にパラメータ指定あり。

元ページは[[fastqc:https://www.bioinformatics.babraham.ac.uk/projects/fastqc/]]

*Trimmomaticによるトリミング [#df89e2ab]
Trimmomaticは[[こちら:http://www.usadellab.org/cms/?page=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のページ:http://www.usadellab.org/cms/?page=trimmomatic]]と[[マニュアル:http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/TrimmomaticManual_V0.32.pdf]]に書かれている。以下、上記の例について説明する(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.fastq|PairedEndでの入力forwardファイル |
|SRR453566_2.fastq|PairedEndでの入力backwardファイル |
|paired_SRR453566}_1.trim.fastq|PairedEndでの出力paired forwardファイル |
|unpaired_SRR453566}_1.trim.fastq|PairedEndでの出力unpaired forwardファイル |
|paired_SRR453566}_2.trim.fastq|PairedEndでの出力paired backwardファイル |
|unpaired_SRR453566}_2.trim.fastq|PairedEndでの出力unpaired backwardファイル |
これ以降は、個別の除去ステップを指定する。
|ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10|ステップ1でIllumina adapterを除去、TruSeq***はアダプターを記述したfastaファイル、2は最大ミスマッチ数、30は回文アライメント時に2つの隣接リードがどれだけ正確にマッチするかを指定、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:https://github.com/timflutre/trimmomatic/blob/master/adapters/TruSeq3-PE-2.fa]]
を使うことができた。

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