![[Toho University]](image/CLogo_Kihon_process.jpg "[Toho University]")
|
Pythonバイオ/ツール/KakenRNAseqhttps://pepper.is.sci.toho-u.ac.jp:443/pepper/index.php?Python%A5%D0%A5%A4%A5%AA%2F%A5%C4%A1%BC%A5%EB%2FKakenRNAseq |
![[YYLAB]](image/yamanouc2.png "[YYLAB]")
|
|
|
Pythonバイオ/ツール/KakenRNAseqhttps://pepper.is.sci.toho-u.ac.jp:443/pepper/index.php?Python%A5%D0%A5%A4%A5%AA%2F%A5%C4%A1%BC%A5%EB%2FKakenRNAseq |
|
参照説明：リボソームRNAの除去
手順
cd /home/yamanouc/src/KishimotoRNA2 for f in 181012_10B_S7_R1_001 181012_1p2-1_S11_R1_001 181012_1p2-2_S13_R1_001 181012_2p5-1_S12_R1_001 181012_2p6-1_S14_R1_001 181012_45a_plus_S2_R1_001 181012_45b_plus_S3_R1_001 181012_45c_plus_S4_R1_001 181013_45A_minus_S18_R1_001 181019_43B_S1_R1_001 181019_45A_plus_S5_R1_001 181019_45L_S6_R1_001 181019_45a10D_plus_S9_R1_001 181019_45aIII6c_plus_S8_R1_001 181019_45d7B_plus_S10_R1_001 181026_10D_minus_S21_R1_001 181026_2-10B_minus_S20_R1_001 181026_45a_minus_S15_R1_001 181026_45b_minus_S16_R1_001 181026_45c_minus_S17_R1_001 181103_45alll6c_minus_S22_R1_001 181103_45d7B_minus_S23_R1_001 181103_Anc_S19_R1_001 PwOw_minus_S25_R1_001 PwOw_plus_S24_R1_001 do nohup sortmerna --ref ./rRNA_databases/silva-bac-16s-id90.fasta,./index/silva-bac-16s-db:\ ./rRNA_databases/silva-bac-23s-id98.fasta,./index/silva-bac-23s-db:\ ./rRNA_databases/rfam-5s-database-id98.fasta,./index/rfam-5s-db:\ ./rRNA_databases/rfam-5.8s-database-id98.fasta,./index/rfam-5.8s-db \ --reads 1st/1st_$f.fastq --sam --num_alignments 1 --fastx --aligned 1st/1st_$f.rRNA \ --other 1st/1st_$f.non_rRNA --log -v > 1st/1st_$f.out & # nohup sortmerna --ref ./rRNA_databases/silva-bac-16s-id90.fasta,./index/silva-bac-16s-db:\ ./rRNA_databases/silva-bac-23s-id98.fasta,./index/silva-bac-23s-db:\ ./rRNA_databases/rfam-5s-database-id98.fasta,./index/rfam-5s-db:\ ./rRNA_databases/rfam-5.8s-database-id98.fasta,./index/rfam-5.8s-db \ --reads 2nd/2nd_$f.fastq --sam --num_alignments 1 --fastx --aligned 2nd/2nd_$f.rRNA \ --other 2nd/2nd_$f.non_rRNA --log -v > 2nd/2nd_$f.out & done
また
cd /home/yamanouc/src/KishimotoRNA3 bash ./sortmerna.sh 3rd_181019_45a10D_plus_S9_R1_001.fastq & bash ./sortmerna.sh 3rd_181019_45d7B_plus_S10_R1_001.fastq & bash ./sortmerna.sh 3rd_181019_45aIII6c_plus_S8_R1_001.fastq & bash ./sortmerna.sh 3rd_181012_2p5-1_S12_R1_001.fastq & bash ./sortmerna.sh 3rd_181012_1p2-1_S11_R1_001.fastq & bash ./sortmerna.sh 3rd_181012_1p2-2_S13_R1_001.fastq & bash ./sortmerna.sh 3rd_181012_45b_plus_S3_R1_001.fastq & bash ./sortmerna.sh 3rd_181019_43B_S1_R1_001.fastq & bash ./sortmerna.sh 3rd_PwOw_plus_S16_R1_001.fastq & bash ./sortmerna.sh 3rd_181012_45c_plus_S4_R1_001.fastq & bash ./sortmerna.sh 3rd_181012_2p6-1_S14_R1_001.fastq & bash ./sortmerna.sh 3rd_181103_Anc_S15_R1_001.fastq & bash ./sortmerna.sh 3rd_181019_45L_S6_R1_001.fastq & bash ./sortmerna.sh 3rd_181019_45A_plus_S5_R1_001.fastq & bash ./sortmerna.sh 3rd_181012_45a_plus_S2_R1_001.fastq &
ただしsortmerna.shは
sortmerna --ref ../KishimotoRNA2/rRNA_databases/silva-bac-16s-id90.fasta,../KishimotoRNA2/index/silva-bac-16s-db:../KishimotoRNA2/rRNA_databases/silva-bac-23s-id98.fasta,../KishimotoRNA2/index/silva-bac-23s-db:../KishimotoRNA2/rRNA_databases/rfam-5s-database-id98.fasta,../KishimotoRNA2/index/rfam-5s-db:../KishimotoRNA2/rRNA_databases/rfam-5.8s-database-id98.fasta,../KishimotoRNA2/index/rfam-5.8s-db --reads $1 --sam --num_alignments 1 --fastx --aligned $1.rRNA --other $1.non_rRNA --log -a 8 -v
結果
readcountのALLのページ
手順
fastqc --nogroup xxx/nonrRNA.fastq
結果
FastQC resultsのページ
それぞれのサンプルのデータは3runとも独立なので、そのままマージしてよいと判断。
具体的には、run-1と2と3をcatで結合した。
⇒出力は<sample>.nonrRNA.fastqとした。
手順
for f in 10B 1p2-1 1p2-2 2p5-1 2p6-1 45a_plus 45b_plus 45c_plus 43B 45A_plus 45L 45a10D_plus 45aIII6c_plus 45d7B_plus Anc PwOw_plus 45A_minus 10D_minus 2-10B_minus 45a_minus 45b_minus 45c_minus 45alll6c_minus 45d7B_minus PwOw_minus do nohup java -jar -Xmx512m /usr/local/Trimmomatic/trimmomatic-0.38.jar SE \ -threads 32 \ -phred33 \ -trimlog $f.trimlog \ $f.nonrRNA.fastq \ $f.nonrRNA.trim.fastq \ ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10 \ LEADING:20 \ TRAILING:20 \ SLIDINGWINDOW:4:15 \ MINLEN:36 > $f.trim.out & done
結果は、
Trimming results のページ
手順は、あらかじめAP012030のインデックスを作っておき、それぞれのサンプルのfastqをマッピング処理する。
Pythonバイオ/ツール/RNAseq�▲蝓璽匹了仮版枸鵑悗離泪奪團鵐� のHISAT2の項参照。
# hisat2.shの中身 for f in 10B 1p2-1 1p2-2 2p5-1 2p6-1 45a_plus 45b_plus 45c_plus 43B 45A_plus 45L 45a10D_plus 45aIII6c_plus 45d7B_plus Anc PwOw_plus 45A_minus 10D_minus 2-10B_minus 45a_minus 45b_minus 45c_minus 45alll6c_minus 45d7B_minus PwOw_minus do nohup hisat2 -p 32 --dta-cufflinks -x AP012030.fna -U $f.nonrRNA.trim.fastq -S $f.sam > $f.hisat2.log & done
結果として、それぞれのSAMファイルができる。hisat2のログ出力<sample>hisat2.logを 集めてリスト化したもの ⇒ hisat2 alignment
その前に、SAMファイルをソートしてかつBAMにしておく。
# compresssam.sh for f in 10B 1p2-1 1p2-2 2p5-1 2p6-1 45a_plus 45b_plus 45c_plus 43B 45A_plus 45L 45a10D_plus 45aIII6c_plus 45d7B_plus Anc PwOw_plus 45A_minus 10D_minus 2-10B_minus 45a_minus 45b_minus 45c_minus 45alll6c_minus 45d7B_minus PwOw_minus do nohup samtools sort -o $f.bam -@ 32 $f.sam > $f.samtools_sort.log & done
得られた <sample>.bam をfeatureCountsでカウントアップする。この時、参照ゲノムの アノテーション GFF ファイルが必要で、かつ、geneフィーチャーに対してgene_idを入れておく必要がある。（geneフィーチャーに入れたことがキーで、後のfeatureCountsの起動パラメータに -t gene を指定する。 -t CDS とすると、CDSフィーチャーにはgene_idが入っていないので、エラーを起こす。
それで作ったアノテーションをAP_e.gffとする。
featureCountsの起動は、
# featureCounts.sh nohup featureCounts -T 32 -t gene -g gene_id -a AP_e.gff -o Counts.txt *.sam > featureCounts.log &
nohupの必要性は無いかも知れない。（それほど時間がかからない）
# geneid_product.ipynb
import pandas as pd
import numpy as np
import urllib.parse
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()
print('N', df.T * 10**3)
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/
print('normalize_tpm')
df_tmp = df.copy()
df_tmp = normalize_per_kilobase(df_tmp, gene_length)
df_tmp = normalize_per_million_reads(df_tmp)
return df_tmp
with open('AP_e.gff', 'r') as fin:
dict = {(0,0): ['','']}
for line in fin:
line=line.rstrip()
if not line.startswith('#'):
s=line.split('\t')
startend = ( int(s[3]), int(s[4]) )
items=s[8].split(';')
for item in items:
if item.startswith('product='):
product = item.split('=')[1]
product = urllib.parse.unquote(product)
if startend not in dict.keys():
dict[startend] = ['', product]
#print('added product -> ', startend, dict[startend])
elif dict[startend][1] != '':
print('dict', startend, 'product part is not empty', dict[startend])
else:
dict[startend] = [dict[startend][0], product]
#print('add product -> ', startend, dict[startend])
if item.startswith('gene_id='):
geneid = item.split('=')[1]
if startend not in dict.keys():
dict[startend] = [geneid, '']
#print('added geneid -> ', startend, dict[startend])
elif dict[startend][0] != '':
print('dict', startend, 'geneid part is not empty', dict[startend])
else:
dict[startend] = [geneid, dict[startend][1]]
#print('add geneid -> ', startend, dict[startend])
fout = open('AP012030_geneid_product.tsv', 'w')
print('%s\t%s' % ('gene_id', 'product'), file=fout)
for k, v in sorted(dict.items(), key=lambda x: x[0][0], reverse=False):
if (v[0]!='' and v[1]!=''):
print('%s\t%s' % (v[0], v[1]), file=fout)
#print(v)
fout.close()
### Counts.txtとマージする
dfc = pd.read_csv('Counts.txt', index_col=0, skiprows=1, sep='\t')
samples = ['10B.sam','10D_minus.sam','1p2-1.sam','1p2-2.sam','2-10B_minus.sam',\
'2p5-1.sam','2p6-1.sam','43B.sam','45A_minus.sam','45A_plus.sam','45L.sam','45a10D_plus.sam',\
'45aIII6c_plus.sam','45a_minus.sam','45a_plus.sam','45alll6c_minus.sam','45b_minus.sam',\
'45b_plus.sam','45c_minus.sam','45c_plus.sam','45d7B_minus.sam','45d7B_plus.sam','Anc.sam',
'PwOw_minus.sam','PwOw_plus.sam']
gene_products = pd.read_csv('AP012030_geneid_product.tsv', index_col=0, sep='\t')
dfg = dfc.join(gene_products, how='inner')
dfg.to_csv('count_prod.csv', sep='\t')
dfs = dfc[samples]
dfs.to_csv('count_raw.csv', sep='\t')
# dfc_fpm = normalize_per_million_reads(dfc)
gene_length = dfc['Length']
#dfs_tpm = normalize_tpm(dfs, gene_length)
dfg_tpm = dfg.copy()
dfg_tpm[samples] = normalize_tpm(dfg[samples], gene_length)
dfg_tpm.to_csv("count_tpm.tsv", sep="\t")
print('complete')
結果のcount_tpm.tsvは このよう になった。
TPM補正後の発現値の分布を確認。
# 比較パターンはC:\Users\yamanouc\Dropbox\RNA重複\KishimotoRNA2\190605 RNAseq サンプル名.xlsx
# Anc - 43B - 45L - *1.2-1 - +2.5-1
# Anc - 43B - 45L - *1.2-2 - +2.6-2
# 45a 2-10B + - PwOw + - 45a III6c + - 45a 10D + - 45d 7B +
# 43B 45a + - 45b + - 45c + 45A + - 45L
# Pythonバイオ/ツール/RNAseq��-2 Python版発現解析（続きの続き）参照
# 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%A4-2%20Python%C8%C7%C8%AF%B8%BD%B2%F2%C0%CF%A1%CA%C2%B3%A4%AD%A4%CE%C2%B3%A4%AD%A1%CB
%matplotlib inline
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()
sns.set_style('whitegrid')
df = pd.read_csv("../count_tpm.tsv", sep="\t", index_col=0)
df = df[['Anc.sam','43B.sam','45L.sam','1p2-1.sam','2p5-1.sam']]
all_zero_index = df.index[df.sum(axis=1) == 0]
df = df.drop(all_zero_index)
dfl = np.log10(df)
dfl = dfl.replace([np.inf, -np.inf], -1) # 値がinf/-infのデータは-1に置換
#dfl.plot.kde(style=['-', '--', ':', '-', '-.', ':']) # 線のパターンは５種類しかない
sns.kdeplot(dfl['Anc.sam']) # 別々に書いて重ねることにした
sns.kdeplot(dfl['43B.sam'])
sns.kdeplot(dfl['45L.sam'])
sns.kdeplot(dfl['1p2-1.sam'])
sns.kdeplot(dfl['2p5-1.sam'])
plt.xlim(-1, 5)
plt.xlabel('$log_{10}(TPM)$')
plt.ylabel('density')
plt.show()
geneごとのサンプル間の変化をヒートマップで表示。全geneを表示するとごちゃごちゃするので、サンプル間の差が小さいgeneは無視することにする。その選択方法は、 変動係数 （標準偏差/平均）が一定値（たとえば0.8）より小さいgeneは無視。
ここでは、0.8を閾値としたところ、267 geneが残った。その残ったgeneをヒートマップ。
%matplotlib inline
# heatmap クラスター化によって寄せる
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.cluster.hierarchy import linkage, dendrogram, leaves_list
from scipy.spatial.distance import pdist
sns.set()
sns.set(style='whitegrid', font=['IPAexGothic'], palette='Spectral')
# count_tpm.tsvの読込み
df = pd.read_csv("../count_tpm.tsv", sep="\t", index_col=0)
df = df[['Anc.sam','43B.sam','45L.sam','1p2-1.sam','2p5-1.sam']]
# サンプルTPM（log前）間の差が小さい行を除く。具体的には変動係数（標準偏差/平均）が0.8以上のものを残す。
dfr = df.copy()
dfr['dev']=(df.std(axis=1)/df.mean(axis=1))
dfr = dfr[dfr['dev'] > 0.8]
print('num dfr', len(dfr))
dfr = dfr + 1 # logを取る前に１を足す
# log10を取る
dfl = np.log10(dfr)
dfl_t = dfl.T
#linkage_result_gene = linkage(dfl, method='average', metric='correlation')
linkage_result_gene = linkage(dfl)
plt.figure(figsize=(16,20))
res = dendrogram(linkage_result_gene, labels=dfl.index, leaf_rotation=90,
leaf_font_size=10, get_leaves=True)
leavesy = leaves_list(linkage_result_gene)
#print('leavesy', leavesy)
lvs = pd.DataFrame([res['ivl'], res['leaves']], columns=res['ivl']).T
lvs.columns = ['gene_id2', 'leaves']
#print(lvs)
dfl2 = pd.merge(lvs, dfl, left_on='gene_id2', right_index=True)
#lvsの並び順のままの方がいい
print(dfl2.head())
plt.figure(figsize=(10, 50))
sns.set_palette('hot')
sns.heatmap(dfl2[['Anc.sam','43B.sam','45L.sam','1p2-1.sam','2p5-1.sam']], vmax=3, vmin=0, center=1.0, cmap='YlOrBr')
plt.savefig('comp1.png')
plt.show()
kdplot.png
