ノート/岸本5/手順メモ～RNAseq

サイトトップ

ノート
データ処理
├マイニング
├バスケット解析
├バスケット解析をＲで
├図書貸出をＲで
├モデルクラスタリング
└
テキストマイニング
├テキストマイニングTM
├TMとシソーラス
├PubMedをTM
├TMとMeSH
├TMとNLTK
├テキストマイニングとtagger
├医薬品添付文書ＤＢ
├論文の処理1
├剽窃１
├Bing検索
└
分子進化学
├分子進化学
└
Pythonと論文アクセス
├Pythonを使ってみる
├PythonでPubMed
├Pythonで...続き
└
Macメモ
├Script

&ref(): File not found: "favicon_story.txt" at page "ノート/ノート";

自信ないので
├Winのバッチ
│├スクリプト入門
│├バッチファイル
│├SETコマンド
│├自分のメモ
│└
└

メニュー編集

岸本　RNAseq　手順をメモっておく 2019-06-04 †

科研の分、ゲノムresequenceデータ

↑

受け取ったデータをコピー †

~/src/KishimotoRNA2/下

↑

結合 †

~/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による品質チェック †

~/src/KishimotoRNA2/下

fastqc --nogroup 181012_10B_S7_R1_001.fastq

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

↑

トリミング †

Trimmomaticでトリミングを行う。ここと同じ処理だが、今回は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

↑

マッピング †

hisat-2を使う。以前にここでテスト済み。

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、準備と実行 †

準備は、gene_idを付けること　ここを参照

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ではないので）付けない。

このあと、メモではpythonでプログラムを書いてカウントを正規化し、適宜処理している。

もしfeatureCountではなくてカウントしたいとき、HISAT-2ではここにメモしてあるように、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
以下略

↑

ポストプロセシング †

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"])
#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 での除算を防ぐため、分&#11935;に微&#12073;な値を加えている
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')

最新の5件

岸本 RNAseq 手順をメモっておく 2019-06-04 †

受け取ったデータをコピー †

結合 †

fastqcによる品質チェック †

トリミング †

マッピング †

FeatureCount、準備と実行 †

ポストプロセシング †

岸本　RNAseq　手順をメモっておく 2019-06-04 †