[[Pythonバイオ/メニュー]]

** ribosomal RNA のフィルタリング [#y5d4c708]
参照説明:[[リボソームRNAの除去>Pythonバイオ/ツール/RNAseq-X リボソーム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:http://pepper.is.sci.toho-u.ac.jp/KakenRNAseq/readcount.htm]]のALLのページ

** FastQCによるQuality Check [#w575adbc]
手順~
 fastqc --nogroup xxx/nonrRNA.fastq

結果~
[[FastQC results:http://pepper.is.sci.toho-u.ac.jp/KakenRNAseq/fastqc/]]のページ


**3つのrunのfastqファイルを結合 [#efa2ce79]
それぞれのサンプルのデータは3runとも独立なので、そのままマージしてよいと判断。

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


**Trimmomaticによるトリミング [#c33dc1eb]
手順~
 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:http://pepper.is.sci.toho-u.ac.jp/KakenRNAseq/trimmomatic.htm]] のページ


** HISAT-2によるAP012030へのマッピング [#s7a5fa72]

手順は、あらかじめ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:http://pepper.is.sci.toho-u.ac.jp/KakenRNAseq/hisat2_alignment.htm]]

** カウントアップその1 featureCountsによるカウントアップ [#eabde175]
その前に、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補正 [#u6f74171]

-productの注記追加
-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は [[このよう:https://pepper.is.sci.toho-u.ac.jp/KakenRNAseq/count_tpm.htm]] になった。

**発現値分布の確認 [#e0d9ac05]
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()

&ref(kdplot.png);

**ヒートマップ [#rec11675]

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()

[[図:]]
[[図:https://pepper.is.sci.toho-u.ac.jp/KakenRNAseq/comp1-08.png]]

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