[[山内のサイト]]
*岸本 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
この２つの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を指定する。入力ファイルは１つだけ、出力ファイルも１つだけ（unpairedは発生しないのでpaired１つだけ）。またリード長がどうも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（Ｒ）で処理する手がある。
後の処理のために、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')
ステップ２
ここから正規化処理。 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それぞれの１、２、３の平均値を求める~
×geneごとに発現変動をlog2 foldとして求める： log2(chemostat/batch)~
×geneごとに３つ並べた表を作る 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')
![[Toho University]](image/CLogo_Kihon_process.jpg "[Toho University]")
![[YYLAB]](image/yamanouc2.png "[YYLAB]")
