star_rsem

**文書の過去の版を表示しています。**

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(最新)が落とせっるぽい。ばんざい

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%

マッピング率がテキストのバージョンよりわずかに向上していた。

RSEMで発現量計算

リファレンス作成

リファレンスを作成。テキスト通り。 

(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!
227301 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!

結果のまとめ

(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

テキストのメインはここまで。あとはIGVでマッピング状況を見てみたり発現量解析に進んだりしましょう。

  • star_rsem.1576050367.txt.gz
  • 最終更新: 2019/12/11 07:46
  • by 133.11.144.10