差分

このページの2つのバージョン間の差分を表示します。

--- star_rsem [2019/12/11 07:46] – [各種ツールインストール] 133.11.144.10
+++ star_rsem [Unknown date] (現在) – 削除 - 外部編集 (Unknown date) 127.0.0.1
@@ 行 1: / 行 1: @@
-## anaconda install
-<code>
-[kijima.yusuke@m48 work]$ pwd
-/home/kijima.yusuke/work
-[kijima.yusuke@m48 work]$ cd work/tool/download/
-[kijima.yusuke@m48 work]$ wget https://repo.anaconda.com/archive/Anaconda3-2019.10-Linux-x86_64.sh
-[kijima.yusuke@m48 work]$ chmod +x Anaconda3-2019.10-Linux-x86_64.sh
-[kijima.yusuke@m48 work]$ ./Anaconda3-2019.10-Linux-x86_64.sh
-[kijima.yusuke@m48 download]$ ./Anaconda3-2019.10-Linux-x86_64.sh
-Welcome to Anaconda3 2019.10
-In order to continue the installation process, please review the license
-agreement.
-Please, press ENTER to continue
->>>
-===================================
-Anaconda End User License Agreement
-===================================
-Copyright 2015, Anaconda, Inc.
-All rights reserved under the 3-clause BSD License:
-...
-Do you accept the license terms? [yes|no]
-[no] >>> yes
-Anaconda3 will now be installed into this location:
-/home/kijima.yusuke/anaconda3
-  - Press ENTER to confirm the location
-  - Press CTRL-C to abort the installation
-  - Or specify a different location below
-[/home/kijima.yusuke/anaconda3] >>> /home/kijima.yusuke/work/tool/anaconda3
-PREFIX=/home/kijima.yusuke/work/tool/anaconda3
-Unpacking payload ...
-...
-Thank you for installing Anaconda3!
-</code>
-一回サーバーから抜ける or bashrcを読み込むことでanaconda環境が完成。ユーザー名の左に(base)の文字が追加される。
-<code>
-(base) [kijima.yusuke@m48 ~]$
-</code>
-## anaconda python3環境の作成
-一応環境を隔離しておく。
-<code>
-(base) [kijima.yusuke@m48 work]$ conda create -n anac_py37 python=3.7 anaconda
-</code>
-終わったら環境を起こす。
-<code>
-(base) [kijima.yusuke@m48 work]$ conda activate anac_py37
-(anac_py37) [kijima.yusuke@m48 work]$ #baseからanac_py37に環境が変わった
-</code>
-## 各種ツールインストール
-テキスト通り。
-<code bash>
-#trimmomatic
-(anac_py37) [kijima.yusuke@m48 work]$ conda install -c bioconda trimmomatic
-#fastqc
-(anac_py37) [kijima.yusuke@m48 work]$ conda install -c bioconda fastqc
-#STAR
-(anac_py37) [kijima.yusuke@m48 work]$ conda install -c bioconda star
-#RSEM
-(anac_py37) [kijima.yusuke@m48 work]$ conda install -c bioconda rsem
-(anac_py37) [kijima.yusuke@m48 work]$ rsem-calculate-expression --version
-Current version: RSEM v1.3.1
-</code>
-テキストとは異なり、現在はconda経由でRSEM v1.3(最新)が落とせっるぽい。ばんざい
-## IGVのインストール
-conda installでIGVの最新版が落とせないか一応調べる。
-<code>
-(anac_py37) [kijima.yusuke@m48 work]$ conda search -c bioconda igv
-Loading channels: done
-# Name                       Version           Build  Channel
-igv                           2.3.98               0  bioconda
-igv                            2.4.6               0  bioconda
-igv                            2.4.9               0  bioconda
-igv                            2.4.9               1  bioconda
-igv                           2.4.16               0  bioconda
-igv                           2.4.17               0  bioconda
-igv                            2.5.2               0  bioconda
-</code>
-最新版の2.6.1はbiocondaにはないっぽい。テキスト通りバイナリを落とす。
-<code>
-(anac_py37) [kijima.yusuke@m48 work]$ mkdir rnaseq-training
-(anac_py37) [kijima.yusuke@m48 work]$ cd rnaseq-training
-(anac_py37) [kijima.yusuke@m48 rnaseq-training]$ mkdir tool
-(anac_py37) [kijima.yusuke@m48 rnaseq-training]$ cd tool
-(anac_py37) [kijima.yusuke@m48 tool]$ wget https://data.broadinstitute.org/igv/projects/downloads/2.7/IGV_Linux_2.7.2.zip
-(anac_py37) [kijima.yusuke@m48 tool]$ unzip IGV_Linux_2.7.2.zip
-</code>
-最新が今(2019/12/9)は2.7.2みたいなのでそれを落とした。
-## リファレンスのダウンロード
-<code>
-(anac_py37) [kijima.yusuke@m48 tool]$ cd ~/work/rnaseq-training
-(anac_py37) [kijima.yusuke@m48 rnaseq-training]$ mkdir human
-(anac_py37) [kijima.yusuke@m48 rnaseq-training]$ mkdir human/ref
-#Genome
-(anac_py37) [kijima.yusuke@m48 rnaseq-training]$ wget -P ~/work/rnaseq-training/human/ref/ ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_32/GRCh38.primary_assembly.genome.fa.gz
-#Annotation
-(anac_py37) [kijima.yusuke@m48 rnaseq-training]$ wget -P ~/work/rnaseq-training/human/ref/ ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_32/gencode.v32.annotation.gtf.gz
-</code>
-## シーケンスデータのダウンロード
-URLのリストを作成。ここはテキストにないので好きにすると良い。逐次ダウンロードでも可。
-<code bash>
-(anac_py37) [kijima.yusuke@m48 rnaseq-training]$ for i in $(seq 1 6); do for j in $(seq 1 2); do echo "ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR718/00${i}/SRR718955${i}/SRR718955${i}_${j}.fastq.gz" ;done; done > ~/work/rnaseq-training/human/data/SRR_Acc_list_forDownload.txt
-(anac_py37) [kijima.yusuke@m48 rnaseq-training]$ cat ~/work/rnaseq-training/human/data/SRR_Acc_list_forDownload.txt
-ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR718/001/SRR7189551/SRR7189551_1.fastq.gz
-ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR718/001/SRR7189551/SRR7189551_2.fastq.gz
-ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR718/002/SRR7189552/SRR7189552_1.fastq.gz
-ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR718/002/SRR7189552/SRR7189552_2.fastq.gz
-ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR718/003/SRR7189553/SRR7189553_1.fastq.gz
-ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR718/003/SRR7189553/SRR7189553_2.fastq.gz
-ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR718/004/SRR7189554/SRR7189554_1.fastq.gz
-ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR718/004/SRR7189554/SRR7189554_2.fastq.gz
-ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR718/005/SRR7189555/SRR7189555_1.fastq.gz
-ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR718/005/SRR7189555/SRR7189555_2.fastq.gz
-ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR718/006/SRR7189556/SRR7189556_1.fastq.gz
-ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR718/006/SRR7189556/SRR7189556_2.fastq.gz
-</code>
-wgetでダウンロード。アクセス切れたら中断されてしまうのでscreen使うなりリモートデスクトップにつなぐなりする。回線に絶対の自信があるなら何も考えなくても良い。それからforループで外部に連続的にアクセスするのはあまりよろしくないのでsleepを挟む。
-<code bash>
-(base) [kijima.yusuke@m48 rnaseq-training]$ for i in $(cat ~/files/m48/kijima.yusuke/work/rnaseq-training/human/data/SRR_Acc_list_forDownload.txt); do wget -P ~/files/m48/kijima.yusuke/work/rnaseq-training/human/data/ $i; sleep 10;  done
-</code>
-後々のためアクセッションIDのリストを作成。
-<code bash>
-(anac_py37) [kijima.yusuke@m48 rnaseq-training]$ mkdir human/data
-(anac_py37) [kijima.yusuke@m48 rnaseq-training]$ for i in SRR718955{1..6}; do echo $i; done > ~/work/rnaseq-training/human/data/SRR_Acc_list.txt
-(anac_py37) [kijima.yusuke@m48 rnaseq-training]$ cat human/data/SRR_Acc_list.txt
-SRR7189551
-SRR7189552
-SRR7189553
-SRR7189554
-SRR7189555
-SRR7189556
-</code>
-##Trimmomaticでアダプター除去+クオリティフィルタリング
-テキスト通りanaconda3/share/trimmomaticからアダプター配列を取ってきたいが、このフォルダがなさそうなので探す。
-<code bash>
-(anac_py37) [kijima.yusuke@m48 rnaseq-training]$ find /home/kijima.yusuke/work/tool/anaconda3/ -name *trimmomatic*
-/home/kijima.yusuke/work/tool/anaconda3/envs/anac_py37/share/trimmomatic
-/home/kijima.yusuke/work/tool/anaconda3/envs/anac_py37/share/trimmomatic-0.39-1
-/home/kijima.yusuke/work/tool/anaconda3/envs/anac_py37/share/trimmomatic-0.39-1/trimmomatic.jar
-/home/kijima.yusuke/work/tool/anaconda3/envs/anac_py37/share/trimmomatic-0.39-1/trimmomatic
-/home/kijima.yusuke/work/tool/anaconda3/envs/anac_py37/conda-meta/trimmomatic-0.39-1.json
-/home/kijima.yusuke/work/tool/anaconda3/envs/anac_py37/bin/trimmomatic
-/home/kijima.yusuke/work/tool/anaconda3/pkgs/trimmomatic-0.39-1.tar.bz2
-/home/kijima.yusuke/work/tool/anaconda3/pkgs/trimmomatic-0.39-1
-/home/kijima.yusuke/work/tool/anaconda3/pkgs/trimmomatic-0.39-1/share/trimmomatic
-/home/kijima.yusuke/work/tool/anaconda3/pkgs/trimmomatic-0.39-1/share/trimmomatic-0.39-1
-/home/kijima.yusuke/work/tool/anaconda3/pkgs/trimmomatic-0.39-1/share/trimmomatic-0.39-1/trimmomatic.jar
-/home/kijima.yusuke/work/tool/anaconda3/pkgs/trimmomatic-0.39-1/share/trimmomatic-0.39-1/trimmomatic
-/home/kijima.yusuke/work/tool/anaconda3/pkgs/trimmomatic-0.39-1/info/recipe/0.32/trimmomatic.sh
-/home/kijima.yusuke/work/tool/anaconda3/pkgs/trimmomatic-0.39-1/info/recipe/0.33/trimmomatic.sh
-/home/kijima.yusuke/work/tool/anaconda3/pkgs/trimmomatic-0.39-1/info/recipe/trimmomatic.py
-/home/kijima.yusuke/work/tool/anaconda3/pkgs/trimmomatic-0.39-1/info/recipe/0.35/trimmomatic.sh
-/home/kijima.yusuke/work/tool/anaconda3/pkgs/trimmomatic-0.39-1/bin/trimmomatic
-</code>
-/home/kijima.yusuke/work/tool/anaconda3/pkgs/trimmomatic-0.39-1/share/の中に入っているようだったので、ここのアダプター配列にリンクを貼る。
-<code>
-(anac_py37) [kijima.yusuke@m48 rnaseq-training]$ ln -s /home/kijima.yusuke/work/tool/anaconda3/pkgs/trimmomatic-0.39-1/share/trimmomatic-0.39-1/adapters/TruSeq3-PE-2.fa human/data/
-(anac_py37) [kijima.yusuke@m48 rnaseq-training]$ ls human/data/
-SRR_Acc_list.txt  SRR_Acc_list_forDownload.txt  TruSeq3-PE-2.fa
-</code>
-Trimmomaticを実行。
-<code bash>
-(anac_py37) [kijima.yusuke@m48 data]$ for id in `cat SRR_Acc_list.txt`;do trimmomatic PE -threads 8 -phred33 -trimlog trimlog_${id}.txt -summary summari_${id}.txt ${id}_1.fastq.gz ${id}_2.fastq.gz ${id}_trim_pair_1.fastq.gz ${id}_trim_unpair_1.fastq.gz ${id}_trim_pair_2.fastq.gz ${id}_trim_unpair_2.fastq.gz ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10 LEADING:20 TRAILING:20 SLIDINGWINDOW:4:15 MINLEN:50; done
-</code>
-終わったらsummrayを見てみる。
-<code bash>
-(base) [kijima.yusuke@m48 data]$ cat summary_SRR7189551.txt
-Input Read Pairs: 30499473
-Both Surviving Reads: 30152544
-Both Surviving Read Percent: 98.86
-Forward Only Surviving Reads: 318460
-Forward Only Surviving Read Percent: 1.04
-Reverse Only Surviving Reads: 25447
-Reverse Only Surviving Read Percent: 0.08
-Dropped Reads: 3022
-Dropped Read Percent: 0.01
-</code>
-テキストと完全一致。すごい。
-##Fastqcでクオリティチェック
-テキスト通り。
-<code bash>
-(anac_py37) [kijima.yusuke@m48 data]$ for id in `cat SRR_Acc_list.txt`; do fastqc -t 12 --nogroup -o fastqc -f fastq ${id}_trim_pair_1.fastq.gz ${id}_trim_pair_2.fastq.gz; done
-</code>
-結果は確認しておきましょう。
-##STARでマッピング
-いよいよマッピング。まずはSTAR用のリファレンス作成から。
-###リファレンス作成
-gunzipファイルを解凍する。テキスト通りgtfは少し編集する。
-<code bash>
-(anac_py37) [kijima.yusuke@m48 data]$ cd ../ref/
-(anac_py37) [kijima.yusuke@m48 ref]$ ls
-GRCh38.primary_assembly.genome.fa.gz  gencode.v32.annotation.gtf.gz
-(anac_py37) [kijima.yusuke@m48 ref]$ gunzip GRCh38.primary_assembly.genome.fa.gz
-(anac_py37) [kijima.yusuke@m48 ref]$ zcat gencode.v32.annotation.gtf.gz | grep -v _PAR_Y > gencode.v32.annotation.gtf
-</code>
-リファレンス作成。オプションの意味はテキスト読んでください。
-<code bash>
-(anac_py37) [kijima.yusuke@m48 ref]$ mkdir star_rsem
-(anac_py37) [kijima.yusuke@m48 ref]$ nohup STAR --runMode genomeGenerate --genomeDir star_rsem --runThreadN 12 --sjdbOverhang 99 --genomeFastaFiles GRCh38.primary_assembly.genome.fa --sjdbGTFfile gencode.v32.annotation.gtf &
-</code>
-###マッピング
-やるだけ
-<code bash>
-(anac_py37) [kijima.yusuke@m48 data]$ for id in `cat SRR_Acc_list.txt`; do mkdir ${id}; STAR --genomeDir ../ref/star_rsem/ --runThreadN 12 --outFileNamePrefix ${id}/ --quantMode TranscriptomeSAM --outSAMtype BAM SortedByCoordinate --readFilesIn ${id}_trim_pair_1.fastq.gz ${id}_trim_pair_2.fastq.gz --readFilesCommand gunzip -c; done
-</code>
-結果を確認。
-<code bash>
-(base) [kijima.yusuke@m48 data]$ cat SRR7189551/Log.final.out
-                          Number of input reads |	30152544
-                      Average input read length |	199
-                                    UNIQUE READS:
-                   Uniquely mapped reads number |	28568001
-                        Uniquely mapped reads % |	94.74%
-                          Average mapped length |	199.28
-                          ...
-                             MULTI-MAPPING READS:
-        Number of reads mapped to multiple loci |	1355491
-             % of reads mapped to multiple loci |	4.50%
-        Number of reads mapped to too many loci |	14845
-             % of reads mapped to too many loci |	0.05%
-</code>
-マッピング率がテキストのバージョンよりわずかに向上していた。
-##RSEMで発現量計算
-###リファレンス作成
-リファレンスを作成。テキスト通り。
-<code bash>
-(base) [kijima.yusuke@m48 data]$ cd ../ref/
-(base) [kijima.yusuke@m48 ref]$ conda activate anac_py37
-(anac_py37) [kijima.yusuke@m48 ref]$ rsem-prepare-reference --gtf gencode.v32.annotation.gtf -p 12 GRCh38.primary_assembly.genome.fa star_rsem/hg
-rsem-extract-reference-transcripts star_rsem/hg 0 gencode.v32.annotation.gtf None 0 GRCh38.primary_assembly.genome.fa
-Parsed 200000 lines
-Parsed 400000 lines
-Parsed 600000 lines
-Parsed 800000 lines
-Parsed 1000000 lines
-Parsed 1200000 lines
-Parsed 1400000 lines
-Parsed 1600000 lines
-Parsed 1800000 lines
-Parsed 2000000 lines
-Parsed 2200000 lines
-Parsed 2400000 lines
-Parsed 2600000 lines
-Parsed 2800000 lines
-Parsing gtf File is done!
-GRCh38.primary_assembly.genome.fa is processed!
-transcripts are extracted.
-Extracting sequences is done!
-Group File is generated!
-Transcript Information File is generated!
-Chromosome List File is generated!
-Extracted Sequences File is generated!
-rsem-preref star_rsem/hg.transcripts.fa 1 star_rsem/hg
-Refs.makeRefs finished!
-Refs.saveRefs finished!
-star_rsem/hg.idx.fa is generated!
-star_rsem/hg.n2g.idx.fa is generated!
-</code>
-###結果のまとめ
-<code bash>
-(anac_py37) [kijima.yusuke@m48 data]$ cd rsem/
-(anac_py37) [kijima.yusuke@m48 rsem]$ rsem-generate-data-matrix SRR718955{1..6}.genes.results | sed -e 's/.genes.results//g' -e 's/"//g' > genes_counts_matrix.tsv
-(anac_py37) [kijima.yusuke@m48 rsem]$ rsem-generate-data-matrix SRR718955{1..6}.isoforms.results | sed -e 's/.isoforms.results//g' -e 's/"//g' > isoforms_counts_matrix.tsv
-(anac_py37) [kijima.yusuke@m48 rsem]$ sed 's/\t/,/g' genes_counts_matrix.tsv > genes_counts_matrix.csv
-(anac_py37) [kijima.yusuke@m48 rsem]$ sed 's/\t/,/g' isoforms_counts_matrix.tsv > isoforms_counts_matrix.csv
-</code>
-テキストのメインはここまで。あとはIGVでマッピング状況を見てみたり発現量解析に進んだりしましょう。