[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
[[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を取り出して、１つずつアクセス
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のページ: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は回文アライメント時に２つの隣接リードがどれだけ正確にマッチするかを指定、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]]
を使うことができた。
終了行:
[[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を取り出して、１つずつアクセス
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のページ: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は回文アライメント時に２つの隣接リードがどれだけ正確にマッチするかを指定、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]]
を使うことができた。
ページ名:
![[Toho University]](image/CLogo_Kihon_process.jpg "[Toho University]")
![[YYLAB]](image/yamanouc2.png "[YYLAB]")
