**文書の過去の版を表示しています。**
anaconda install
[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!
一回サーバーから抜ける or bashrcを読み込むことでanaconda環境が完成。ユーザー名の左に(base)の文字が追加される。
(base) [kijima.yusuke@m48 ~]$
anaconda python3環境の作成
一応環境を隔離しておく。
(base) [kijima.yusuke@m48 work]$ conda create -n anac_py37 python=3.7 anaconda
終わったら環境を起こす。
(base) [kijima.yusuke@m48 work]$ conda activate anac_py37 (anac_py37) [kijima.yusuke@m48 work]$ #baseからanac_py37に環境が変わった
各種ツールインストール
テキスト通り。
#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
テキストとは異なり、現在はconda経由でRSEM v1.3(最新)が落とせっるぽい。ばんざい </code>
IGVのインストール
conda installでIGVの最新版が落とせないか一応調べる。
(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
最新版の2.6.1はbiocondaにはないっぽい。テキスト通りバイナリを落とす。
(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
最新が今(2019/12/9)は2.7.2みたいなのでそれを落とした。
リファレンスのダウンロード
(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
シーケンスデータのダウンロード
URLのリストを作成。ここはテキストにないので好きにすると良い。逐次ダウンロードでも可。
(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
wgetでダウンロード。アクセス切れたら中断されてしまうのでscreen使うなりリモートデスクトップにつなぐなりする。回線に絶対の自信があるなら何も考えなくても良い。それからforループで外部に連続的にアクセスするのはあまりよろしくないのでsleepを挟む。
(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
後々のためアクセッションIDのリストを作成。
(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
Trimmomaticでアダプター除去+クオリティフィルタリング
テキスト通りanaconda3/share/trimmomaticからアダプター配列を取ってきたいが、このフォルダがなさそうなので探す。
(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
/home/kijima.yusuke/work/tool/anaconda3/pkgs/trimmomatic-0.39-1/share/の中に入っているようだったので、ここのアダプター配列にリンクを貼る。
(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
Trimmomaticを実行。
(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
終わったらsummrayを見てみる。
(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
テキストと完全一致。すごい。
Fastqcでクオリティチェック
テキスト通り。
(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
結果は確認しておきましょう。
STARでマッピング
いよいよマッピング。まずはSTAR用のリファレンス作成から。
リファレンス作成
gunzipファイルを解凍する。テキスト通りgtfは少し編集する。
(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
リファレンス作成。オプションの意味はテキスト読んでください。
(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 &
マッピング
やるだけ
(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
結果を確認。
(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%
マッピング率がテキストのバージョンよりわずかに向上していた。