[[山内のサイト]]

*岸本 RNAseq 手順をメモっておく 2019-06-04 [#gf7535db]
科研の分、ゲノムresequenceデータ


**受け取ったデータをコピー [#g0fde037]
~~/src/KishimotoRNA2/下

**結合 [#z7c9facf]
~~/src/KishimotoRNA2/下

1st下の1st/1st_181012_10B_S7_R1_001.fastq.gzを解凍し
 gunzip: 1st/1st_181012_10B_S7_R1_001.fastq.gz
2nd下の2nd/2nd_181012_10B_S7_R1_001.fastq.gzを解凍し
 gunzip: 2nd/2nd_181012_10B_S7_R1_001.fastq.gz
この2つのfastqを結合する
 cat: 1st/1st_181012_10B_S7_R1_001.fastq 2nd/2nd_181012_10B_S7_R1_001.fastq > 181012_10B_S7_R1_001.fastq

(これを含めて、以下の処理を、すべてのサンプルについて行う)

**fastqcによる品質チェック [#kc1a5cab]
~~/src/KishimotoRNA2/下

 fastqc --nogroup 181012_10B_S7_R1_001.fastq

⇒181012_10B_S7_R1_001_fastq.html と181012_10B_S7_R1_001_fastqc.zip生成


**トリミング [#qa2c0468]
Trimmomaticでトリミングを行う。[[ここ:https://pepper.is.sci.toho-u.ac.jp/pepper/index.php?Python%A5%D0%A5%A4%A5%AA%2F%A5%C4%A1%BC%A5%EB%2FRNAseq%AD%A1%A5%EA%A1%BC%A5%C9%A5%C7%A1%BC%A5%BF%BC%E8%C6%C0%A4%C8%A5%AF%A5%AA%A5%EA%A5%C6%A5%A3%A5%C1%A5%A7%A5%C3%A5%AF%A1%A6%A5%C8%A5%EA%A5%DF%A5%F3%A5%B0#df89e2ab]]と同じ処理だが、今回はPair EndではなくてSingle Endなので、指定が違う。SEを指定する。入力ファイルは1つだけ、出力ファイルも1つだけ(unpairedは発生しないのでpaired1つだけ)。またリード長がどうも50なので、MINLENは36ぐらいがよかろう。

 java -jar -Xmx4g /usr/local/Trimmomatic/trimmomatic-0.38.jar SE -threads 16 -phred33 181012_10B_S7_R1_001.fastq  181012_10B_S7_R1_001.trim.fastq -trimlog trimlog_181012_10B_S7_R1_001.log ILLUMINACLIP:/usr/local/Trimmomatic/adapters/TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

**マッピング [#y90fbf5b]

hisat-2を使う。[[以前にここ>Pythonバイオ/ツール/RNAseq〜AP012030]]でテスト済み。

 hisat2-build AP012030-new.fasta AP012030.fna (多分)

出力は AP012030.fna.*.ht2 なのでこれらのファイルをコピーして持ってくる。
 AP012030.fna
 AP012030.fna.1.ht2
 AP012030.fna.2.ht2
 AP012030.fna.3.ht2
 AP012030.fna.4.ht2
 AP012030.fna.5.ht2
 AP012030.fna.6.ht2
 AP012030.fna.7.ht2
 AP012030.fna.8.ht2

 hisat2 -p32 -x AP012030.fna -U 181012_10B_S7_R1_001.trim.fastq -S 181012_10B_S7_R1_001.sam

**FeatureCount、準備と実行 [#q5a75af0]
準備は、gene_idを付けること [[ここ>https://pepper.is.sci.toho-u.ac.jp/pepper/index.php?Python%A5%D0%A5%A4%A5%AA%2F%A5%C4%A1%BC%A5%EB%2FRNAseq%AD%A2%A5%EA%A1%BC%A5%C9%A4%CE%BB%B2%BE%C8%C7%DB%CE%F3%A4%D8%A4%CE%A5%DE%A5%C3%A5%D4%A5%F3%A5%B0#qaa4366c]] を参照

 pip install bcbio-gff
がしてある環境で
 python /usr/local/RNAsq/add_gene_id2.py AP012030.gff AP012030_e.gff
但し、
 ~/src/RNAseq-Saccha/遺伝研/2-1_and_2-2/add_gene_id
を改造したバージョン
 /usr/local/RNAseq/add_gene_id2.py
を使う。gene_idの番号の後ろにgene名を追加しかつCDSにもgene_id番号を付ける。

また、前ステージで出来たSAMファイルをBAMファイルにし、次にソートする。
 samtools view -@ 32 -b SRR453566.sam -o SRR453566.bam
 samtools sort -@ 32 SRR453566.bam > SRR453566.sorted.bam 

これに対して、[[featureCounthttps://bi.biopapyrus.jp/rnaseq/expression/featurecounts.html]]を使って
 featureCounts -T 8 -t exon -g gene_id -a AP012030.gff -o  AP012030.counts.txt  *.sorted.bam
このとき、-pは(paired endではないので)付けない。

このあと、[[メモ:https://pepper.is.sci.toho-u.ac.jp/pepper/index.php?Python%A5%D0%A5%A4%A5%AA%2F%A5%C4%A1%BC%A5%EB%2FRNAseq%AD%A3-0%A1%CA%B0%E4%C5%C1%B8%A6%A5%BB%A5%DF%A5%CA%A1%BC%A4%C7%A4%CE%A1%CB%C8%AF%B8%BD%B2%F2%C0%CF%A1%C1%A5%AB%A5%A6%A5%F3%A5%C8%A5%C7%A1%BC%A5%BF%A4%CE%B2%F2%C0%CF#m0a943f2]]ではpythonでプログラムを書いてカウントを正規化し、適宜処理している。


もしfeatureCountではなくてカウントしたいとき、HISAT-2では[[ここ:https://pepper.is.sci.toho-u.ac.jp/pepper/index.php?Python%A5%D0%A5%A4%A5%AA%2F%A5%C4%A1%BC%A5%EB%2FRNAseq%AD%A2%A5%EA%A1%BC%A5%C9%A4%CE%BB%B2%BE%C8%C7%DB%CE%F3%A4%D8%A4%CE%A5%DE%A5%C3%A5%D4%A5%F3%A5%B0#e7486864]]にメモしてあるように、StringTie⇒Ballgawn(R)で処理する手がある。


後の処理のために、gene_idとproductの対応表 gene_id_product.csv を作っておく。
 import pandas as pd
 import numpy as np
 import urllib.parse
 outputlist = []
 with open('AP012030_e.gff', 'r') as fin:
     for line in fin:
         product = ''
         geneid = ''
         line=line.rstrip()
         if not line.startswith('#'):
             s=line.split('\t')
             items=s[8].split(';')
             for item in items:
                 if item.startswith('product='):
                     product = item.split('=')[1]
                     product = urllib.parse.unquote(product)
                 if item.startswith('gene_id='):
                     geneid = item.split('=')[1]
                 if (geneid != '') and (product != ''):
                     outputlist.append(geneid+'\t'+product)
 outputlist = sorted(list(set(outputlist)))
 fout = open('AP012030_gene_id_product.tsv', 'w')
 for u in outputlist:
     print(u, file=fout)
 fout.close()
 print('complete')

表は
 gene_0001_thrL  thr operon leader peptide
 gene_0003_thrA  bifunctional aspartokinase I/homeserine dehydrogenase I
 gene_0005_thrB  homoserine kinase
 gene_0007_thrC  threonine synthase
 gene_0008_yaaX  predicted protein
 gene_0009_yaaA  hypothetical protein
 gene_0011_yaaJ  putative transporter
 gene_0014_talB  transaldolase 1
 以下略

**ポストプロセシング [#sce7fbc5]
featureCountの出力AP012030.counts.txtと、AP012030_gene_id_product.tsvとをマージすることで、AP012030.counts.txtにproductの入った表を作り、AP012030_count_preprocessed.tsvに書き出す。

 import pandas as pd
 import numpy as np
 df = pd.read_table('AP012030.counts.txt', index_col=0, skiprows=1) 
 
 rename_columns = {
 '181012_10B_S7_R1_001.sorted.bam': '45a 2-10B +',
 '181012_1p2-1_S11_R1_001.sorted.bam': '(+1.2)-1',
 '181012_1p2-2_S13_R1_001.sorted.bam': '(+1.2)-2',
 '181012_2p5-1_S12_R1_001.sorted.bam': '(+2.5)-1',
 '181012_2p6-1_S14_R1_001.sorted.bam': '(+2.6)-2',
 '181012_45a_plus_S2_R1_001.sorted.bam': '45a +',
 '181012_45b_plus_S3_R1_001.sorted.bam': '45b +',
 '181012_45c_plus_S4_R1_001.sorted.bam': '45c +',
 '181013_45A_minus_S18_R1_001.sorted.bam': '45A -',
 '181019_43B_S1_R1_001.sorted.bam': '43B',
 '181019_45A_plus_S5_R1_001.sorted.bam': '45A +',
 '181019_45L_S6_R1_001.sorted.bam': '45L',
 '181019_45a10D_plus_S9_R1_001.sorted.bam': '45a 10D +',
 '181019_45aIII6c_plus_S8_R1_001.sorted.bam': '45a III6c +',
 '181019_45d7B_plus_S10_R1_001.sorted.bam': '45d 7B +',
 '181026_10D_minus_S21_R1_001.sorted.bam': '45a 10D -',
 '181026_2-10B_minus_S20_R1_001.sorted.bam': '45a 2-10B -',
 '181026_45a_minus_S15_R1_001.sorted.bam': '45a -',
 '181026_45b_minus_S16_R1_001.sorted.bam': '45b -',
 '181026_45c_minus_S17_R1_001.sorted.bam': '45c -',
 '181103_45alll6c_minus_S22_R1_001.sorted.bam': '45a III6c -',
 '181103_45d7B_minus_S23_R1_001.sorted.bam': '45d 7B -',
 '181103_Anc_S19_R1_001.sorted.bam': 'Anc',
 'PwOw_minus_S25_R1_001.sorted.bam': 'PwOw -',
 'PwOw_plus_S24_R1_001.sorted.bam': 'PwOw +'
 }
 df = df.rename_axis(mapper=rename_columns, axis=1)
 #print(df.head())
 
 # df = df[df.Chr != 'NC_001224.1']  # ミトコンドリア(NC_001224.1)上の遺伝子を除く
 gene_products = pd.read_table('AP012030_gene_id_product.tsv', index_col=0, 
names=["gene_id", "product"])
 gene_products = pd.read_table('AP012030_gene_id_product.tsv', index_col=0, names=["gene_id", "product"])
 #print(gene_products.head())
 
 # マージ
 #df_with_product = gene_products.join(df)  # 値がfloatになる?
 df_with_product = df.join(gene_products)
 #print(df_with_product.head())
 df_count = df_with_product[[
 '45a 2-10B +', '(+1.2)-1', '(+1.2)-2', '(+2.5)-1', '(+2.6)-2', '45a +', '45b +',
 '45c +', '45A -', '43B', '45A +', '45L', '45a 10D +', '45a III6c +', '45d 7B +',
 '45a 10D -', '45a 2-10B -', '45a -', '45b -', '45c -', '45a III6c -', '45d 7B -',
 'Anc', 'PwOw -', 'PwOw +']]
 print(df_count.head())
 df_count.to_csv('AP012030_count_raw.tsv', sep='\t')
 df_with_product.to_csv('AP012030_count_preprocessed.tsv', sep='\t')
 print('complete')



ステップ2

ここから正規化処理。 RPM/FPM、RPKM/FPKM、TPMを計算する。

RPM/FPM Reads per Million / Fragments per Million~
RPKM/FPKM Reads per kilobase of exon per million reads mapped / Fragments per ...~
TPM Transcripts per kilobase million~
FPKMを使って、batch cultureとchemostatの間で発現変動の大きい遺伝子を抽出~

×geneごとにbatch, chemostatそれぞれの1、2、3の平均値を求める~
×geneごとに発現変動をlog2 foldとして求める: log2(chemostat/batch)~
×geneごとに3つ並べた表を作る  batch, chemostat, log2fold。 更にgene_productsを並べる~
×カウントが0のデータを除き、los2fold値でソートする。出力→diff_ex.tsv~

 %matplotlib inline
 # 正規化処理。 RPM/FPM、RPKM/FPKM、TPMを計算する
 # × geneごとに発現変動をlog2 foldとして求める: log2(chemostat/batch)
 # × カウントが0のデータを除き、los2fold値でソートする。出力→diff_ex.tsv
 import pandas as pd
 import numpy as np
 
 def normalize_per_million_reads(df):
     # RPM/FPM = reads per million  fragments per million
     sum_count = df.sum()
     return 10**6 * df / sum_count
 
 def normalize_per_kilobase(df, gene_length):
     # FPKM = fragments per kilobase of exon per million reads mapped
     df_tmp = df.copy()
     df_tmp = (df.T * 10**3 / gene_length).T 
     return df_tmp
 
 def normalize_tpm(df, gene_length):
     # TPM = transcripts per kilobase million   https://bi.biopapyrus.jp/ 
     df_tmp = df.copy()
     df_tmp = normalize_per_kilobase(df_tmp, gene_length)
     df_tmp = normalize_per_million_reads(df_tmp)
     return df_tmp
 
 df_count = pd.read_table('AP012030_count_raw.tsv', index_col=0)
 df_with_product = pd.read_table('AP012030_count_preprocessed.tsv', index_col=0)
 gene_products = pd.read_table('AP012030_gene_id_product.tsv', index_col=0, 
 names=["gene_id", "product"])
 
 # RPM/FPM = reads per million  fragments per million
 df_count_fpm = normalize_per_million_reads(df_count)
 #print(df_count_fpm.head())
 df_count_fpm.to_csv("AP012030_count_fpm.tsv", sep="\t")
 
 # RPKM/FPKM
 # FPKM = fragments per kilobase of exon per million reads mapped
 # single-endの場合、RPKM = reads per kilobase of exon per million reads mapped
 gene_length = df_with_product['Length']
 df_count_fpkm = normalize_per_kilobase(df_count_fpm, gene_length)
 #print(df_count_fpkm.head())
 df_count_fpkm.to_csv('AP012030_count_fpkm.tsv', sep='\t')
 
 # TPM normalization
 # TPM = transcripts per kilobase million   https://bi.biopapyrus.jp/ 
 df_count_tpm = normalize_tpm(df_count, gene_length)
 #print(df_count_tpm.head())
 #print(df_count_tpm.sum())
 df_count_tpm.to_csv("AP012030_count_tpm.tsv", sep="\t")
 
 '''
 # ここからは、今はやらない
 # 発現変動の大きい遺伝子の抽出
 # batch cultureの平均を求める
 df_count_fpkm["batch"] = (df_count_fpkm["batch_1"] + df_count_fpkm["batch_2"] + 
 df_count_fpkm["batch_3"]) / 3
 # chemostatの平均を求める
 df_count_fpkm["chemostat"] = (df_count_fpkm["chemostat_1"] + 
 df_count_fpkm["chemostat_2"] + df_count_fpkm["chemostat_3"]) / 3
 # 発現変動をlog2 fold として求める
 # 0 での除算を防ぐため、分⺟に微⼩な値を加えている
 df_count_fpkm["log2fold"] = \
    df_count_fpkm["chemostat"] / (df_count_fpkm["batch"] + 10**-6)
 df_count_fpkm["log2fold"] = df_count_fpkm["log2fold"].apply(np.log2)
 #print(df_count_fpkm.head())
 # 必要部分のみ抜き出し、productと結合
 diff_ex = df_count_fpkm[["batch", "chemostat", "log2fold"]]
 diff_ex = diff_ex.join(gene_products)
 # カウント数が0のデータを除く
 diff_ex = diff_ex[diff_ex["batch"] > 0]
 diff_ex = diff_ex[diff_ex["chemostat"] > 0]
 #print(diff_ex)
 diff_ex = diff_ex.sort_values("log2fold", ascending=False)
 #print(diff_ex.head())
 #print(diff_ex.tail())
 diff_ex.to_csv("diff_ex.tsv", sep="\t")
 '''
 print('complete')

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