Pythonバイオ/メニュー

ribosomal RNA のフィルタリング

参照説明:リボソーム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によるQuality Check

手順

fastqc --nogroup xxx/nonrRNA.fastq

結果
FastQC resultsのページ

3つのrunのfastqファイルを結合

それぞれのサンプルのデータは3runとも独立なので、そのままマージしてよいと判断。

具体的には、run-1と2と3をcatで結合した。
  ⇒出力は<sample>.nonrRNA.fastqとした。

Trimmomaticによるトリミング

手順

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 のページ

HISAT-2によるAP012030へのマッピング

手順は、あらかじめ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

カウントアップその1 featureCountsによるカウントアップ

その前に、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の必要性は無いかも知れない。(それほど時間がかからない)

featureCountsの結果のTPM補正

# 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=['-', '--', ':', '-', '-.', ':'])  # 線のパターンは5種類しかない
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()

kdplot.png

ヒートマップ

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を取る前に1を足す
# 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()


添付ファイル: filekdplot.png 183件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2019-09-01 (日) 17:00:28 (394d)