################################################ ### Step1-1 ### SRR616268sub_1.fastq.gzの前処理 ### 参考: 前処理 | トリミング | 指定した末端塩基数だけ除去 | 例題4 ### 目的: 3'側のあやしいものを除去 ### 挙動: 1,000,000 reads --> 1,000,000 reads ################################################ in_f <- "SRR616268sub_1.fastq.gz" #入力ファイル名を指定してin_fに格納 out_f <- "hoge1.fastq.gz" #出力ファイル名を指定してout_fに格納 param_trim <- 7 #3'末端のトリムしたい塩基数を指定 #必要なパッケージをロード library(Biostrings) #パッケージの読み込み library(ShortRead) #パッケージの読み込み #入力ファイルの読み込み fastq <- readFastq(in_f) #in_fで指定したファイルの読み込み sread(fastq) #配列情報を表示 #本番 hoge <- width(sread(fastq)) - param_trim#トリム後のend位置情報を取得 hoge[hoge < 1] <- 1 #トリム後のend位置が1未満の場合には1にしている hoge1 <- DNAStringSet(sread(fastq), start=1, end=hoge)#sread(fastq)から指定した範囲を抽出した結果をDNAStringSet形式で格納 hoge2 <- BStringSet(quality(quality(fastq)), start=1, end=hoge)#quality(quality(fastq))から指定した範囲を抽出した結果をBStringSet形式で格納 fastq <- ShortReadQ(hoge1, hoge2, id(fastq))#ShortReadQというクラスオブジェクトを作成してfastqに格納 sread(fastq) #配列情報を表示 quality(fastq) #quality情報を表示 #ファイルに保存 writeFastq(fastq, out_f, compress=T) #fastqの中身を指定したファイル名で保存 ################################################ ### Step1-2 ### SRR616268sub_1.fastq.gzの前処理 ### 参考: 前処理 | トリミング | アダプター配列除去(基礎) | QuasR(Gaidatzis_2015) | 例題5 ### 目的: 5'側の「TruSeq Adapter, Index 3」を除去 ### 挙動: 1,000,000 reads --> 998,659 reads ################################################ in_f <- "hoge1.fastq.gz" #入力ファイル名を指定してin_fに格納 out_f <- "hoge2.fastq.gz" #出力ファイル名を指定してout_fに格納 param_adapter <- "GATCGGAAGAGCACACGTCTGAACTCCAGTCACTTAGGCATCTCGTATGCCGTCTTCTGCTTG"#アダプター配列を指定 param_nrec <- 500000 #一度に処理するリード数を指定 #必要なパッケージをロード library(QuasR) #パッケージの読み込み #本番(アダプター配列除去) res <- preprocessReads(filename=in_f, #in_fを入力とする outputFilename=out_f, #out_fを出力とする Lpattern=param_adapter, #リードの左側(Left)のアダプターを除去 nrec=param_nrec) #一度に処理するリード数 res #確認してるだけです file.size(in_f) #入力ファイルのサイズを表示 file.size(out_f) #出力ファイルのサイズを表示 ################################################ ### Step1-3 ### SRR616268sub_1.fastq.gzの前処理 ### 参考: 前処理 | トリミング | アダプター配列除去(基礎) | QuasR(Gaidatzis_2015) | 例題7 ### 目的: 5'側の「TruSeq Adapter, Index 2」を除去 ### 挙動: 998,659 reads --> 998,649 reads ################################################ in_f <- "hoge2.fastq.gz" #入力ファイル名を指定してin_fに格納 out_f <- "hoge3.fastq.gz" #出力ファイル名を指定してout_fに格納 param_adapter <- "GATCGGAAGAGCACACGTCTGAACTCCAGTCACCGATGTATCTCGTATGCCGTCTTCTGCTTG"#アダプター配列を指定 param_nrec <- 500000 #一度に処理するリード数を指定 #必要なパッケージをロード library(QuasR) #パッケージの読み込み #本番(アダプター配列除去) res <- preprocessReads(filename=in_f, #in_fを入力とする outputFilename=out_f, #out_fを出力とする Lpattern=param_adapter, #リードの左側(Left)のアダプターを除去 nrec=param_nrec) #一度に処理するリード数 res #確認してるだけです file.size(in_f) #入力ファイルのサイズを表示 file.size(out_f) #出力ファイルのサイズを表示 ################################################ ### Step2-1 ### SRR616268sub_2.fastq.gzの前処理 ### 参考: 前処理 | トリミング | 指定した末端塩基数だけ除去 | 例題5 ### 目的: 3'側のあやしいものを除去 ### 挙動: 1,000,000 reads --> 1,000,000 reads ################################################ in_f <- "SRR616268sub_2.fastq.gz" #入力ファイル名を指定してin_fに格納 out_f <- "hoge4.fastq.gz" #出力ファイル名を指定してout_fに格納 param_trim <- 2 #3'末端のトリムしたい塩基数を指定 #必要なパッケージをロード library(Biostrings) #パッケージの読み込み library(ShortRead) #パッケージの読み込み #入力ファイルの読み込み fastq <- readFastq(in_f) #in_fで指定したファイルの読み込み sread(fastq) #配列情報を表示 #本番 hoge <- width(sread(fastq)) - param_trim#トリム後のend位置情報を取得 hoge[hoge < 1] <- 1 #トリム後のend位置が1未満の場合には1にしている hoge1 <- DNAStringSet(sread(fastq), start=1, end=hoge)#sread(fastq)から指定した範囲を抽出した結果をDNAStringSet形式で格納 hoge2 <- BStringSet(quality(quality(fastq)), start=1, end=hoge)#quality(quality(fastq))から指定した範囲を抽出した結果をBStringSet形式で格納 fastq <- ShortReadQ(hoge1, hoge2, id(fastq))#ShortReadQというクラスオブジェクトを作成してfastqに格納 sread(fastq) #配列情報を表示 quality(fastq) #quality情報を表示 #ファイルに保存 writeFastq(fastq, out_f, compress=T) #fastqの中身を指定したファイル名で保存 ################################################ ### Step2-2 ### SRR616268sub_2.fastq.gzの前処理 ### 参考: 前処理 | トリミング | アダプター配列除去(基礎) | QuasR(Gaidatzis_2015) | 例題5 ### 目的: 5'側の「Illumina Single End PCR Primer 1」を除去 ### 挙動: 1,000,000 reads --> 999,233 reads ################################################ in_f <- "hoge4.fastq.gz" #入力ファイル名を指定してin_fに格納 out_f <- "hoge5.fastq.gz" #出力ファイル名を指定してout_fに格納 param_adapter <- "AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT"#アダプター配列を指定 param_nrec <- 500000 #一度に処理するリード数を指定 #必要なパッケージをロード library(QuasR) #パッケージの読み込み library(Biostrings) #パッケージの読み込み #前処理(reverse complementの作成) hoge <- reverseComplement(DNAString(param_adapter))#逆相補鎖(reverse complement)を作成 #本番(アダプター配列除去) res <- preprocessReads(filename=in_f, #in_fを入力とする outputFilename=out_f, #out_fを出力とする Lpattern=as.character(hoge),#リードの左側(Left)のアダプターを除去 nrec=param_nrec) #一度に処理するリード数 res #確認してるだけです file.size(in_f) #入力ファイルのサイズを表示 file.size(out_f) #出力ファイルのサイズを表示 ################################################ ### Step3-1 ### paired-endに戻す ### 参考: 前処理 | フィルタリング | paired-end | 共通リード抽出 | ShortRead(Morgan_2009) ### 目的: リード数の異なる2つのファイルを入力として、共通リードからなるファイルを得る ### 挙動: 998,649 reads and 999,233 reads --> 998,521 reads ################################################ in_f1 <- "hoge3.fastq.gz" #入力ファイル名を指定してin_f1に格納 in_f2 <- "hoge5.fastq.gz" #入力ファイル名を指定してin_f2に格納 out_f1 <- "SRR616268sub_trim3_1.fastq.gz"#出力ファイル名を指定してout_f1に格納(59,092,219 bytes) out_f2 <- "SRR616268sub_trim3_2.fastq.gz"#出力ファイル名を指定してout_f2に格納(54,667,920 bytes) #必要なパッケージをロード library(ShortRead) #パッケージの読み込み #入力ファイルの読み込み fastq1 <- readFastq(in_f1) #in_f1で指定したファイルの読み込み fastq2 <- readFastq(in_f2) #in_f2で指定したファイルの読み込み head(id(fastq1), n=4) #description部分を表示 head(id(fastq2), n=4) #description部分を表示 #前処理(FASTQのdescription行を整形; in_f1) fastq <- fastq1 #fastqとして取り扱う hoge <- strsplit(as.character(id(fastq)), " ", fixed=TRUE)#id(fastq)中の文字列を" "で区切った結果をリスト形式でhogeに格納 description <- BStringSet(sapply(hoge,"[[", 1))#hogeのリスト中の1番目の要素のみ取り出してdescriptionに格納 fastq <- ShortReadQ(sread(fastq), quality(fastq), description)#ReadFastQ関数を用いてReadFastQというクラスオブジェクトを一から作成 head(id(fastq), n=4) #description部分を表示 fastq1 <- fastq #元に戻す #前処理(FASTQのdescription行を整形; in_f2) fastq <- fastq2 #fastqとして取り扱う hoge <- strsplit(as.character(id(fastq)), " ", fixed=TRUE)#id(fastq)中の文字列を" "で区切った結果をリスト形式でhogeに格納 description <- BStringSet(sapply(hoge,"[[", 1))#hogeのリスト中の1番目の要素のみ取り出してdescriptionに格納 fastq <- ShortReadQ(sread(fastq), quality(fastq), description)#ReadFastQ関数を用いてReadFastQというクラスオブジェクトを一から作成 head(id(fastq), n=4) #description部分を表示 fastq2 <- fastq #元に戻す #本番 common <- intersect(id(fastq1), id(fastq2))#共通して現れるリードID情報を抽出 length(fastq1) #fastq1のリード数を表示 length(fastq2) #fastq2のリード数を表示 length(common) #共通リード数を表示 obj1 <- is.element(as.character(id(fastq1)), as.character(common))#共通リードの位置をTRUE、それ以外をFALSEとしたベクトルを作成 obj2 <- is.element(as.character(id(fastq2)), as.character(common))#共通リードの位置をTRUE、それ以外をFALSEとしたベクトルを作成 #ファイルに保存 writeFastq(fastq1[obj1], out_f1, compress=T)#fastqの中身を指定したファイル名で保存 writeFastq(fastq2[obj2], out_f2, compress=T)#fastqの中身を指定したファイル名で保存 file.size(in_f1) #入力ファイルのサイズを表示 file.size(in_f2) #入力ファイルのサイズを表示 file.size(out_f1) #出力ファイルのサイズを表示 file.size(out_f2) #出力ファイルのサイズを表示