NGS alignment

---

1 aligner

---

1.1 hisat2

• アラインメントの結果は-S(sam形式のoutput)を指定しなければ, 標準出力にsam形式として出力される. 標準エラー出力にマップ率の情報が出力される.　samファイルを扱うことはほとんどないのでリダイレクトでsamtoolsに渡してbam形式のファイルを作成.
• -U: シングルエンドの場合fastqファイルの指定
• -1, -2 : ペアエンドの場合のfastqファイルの指定
• -rna-strandness: strand指定 FR, RF,F,R
• --dta: トランスクリプト用に調整されたアライメントの実行(stringtieを実行する場合推奨)
• --no-spliced-alignment:
• --un-conc-gz: ペアが適切にマップされなかったリードのfastqファイル
• --summary-file: print alignment summary to this file.
• --new-summary : print alignment summary in a new style, which is more machine-friendly.

# リファレンス作成
hisat2-build -p 4 ref.fna ref 

# シングル
hisat2 -p 4 -x ref -U seq.fastq.gz > seq.sam
hisat2 -p 4 -x ref -U seq.fastq.gz 2> log.txt | samtools sort -@ 4 -O BAM -o out.sort.bam -

# ペアエンド
hisat2 -p 4 -x ref -1 seq_R1.fastq.gz -2 seq_R2.fastq.gz | samtools sort -@ 4 -O BAM -o out.sort.bam -

# 引数が長くなる場合は改行を入れる
hisat2 -p 4 -x ~/db/index/hisat2_idx/ref \
-1 seq_R1.fastq.gz -2 seq_R2.fq.gz \
--rna-strandness RF \
--un-conc-gz fail.fastq.gz \
--met-file met.txt \
--no-spliced-alignment 2> log.txt | \
samtools sort -@ 4 -O BAM -o out.sort.bam -

1.2 bowtie2

# リファレンス作成 
bowtie2-build -f ref.fa ref

# シングルエンド
bowtie2 -p 4 -x ref seq.fastq.gz

# ペアエンド
bowtie2 -p 4 -x ref -1 R1.fastq.gz -2 R2.fastq.gz

1.3 bwa

• -t: スレッド数 (デファルト1)
• -R: リードのヘッダー情報。
• -M: mark shorter split hits as secondary

# index作成 
bwa index -p ref ref.fa

# bwa mem
bwa mem -t 4 ref seq_R1.fastq.gz | samtools sort -@8 -O BAM -o out.sort.bam -
bwa mem -t 4 ref seq_R1.fastq.gz seq_R2.fastq.gz| samtools sort -@8 -O BAM -o out.sort.bam -

1.4 複数サンプルをループ処理.

• 出力ファイル名などを指定する場合入力ファイル名から部分文字列を取り出す.
  • ubuntuではpfx=${r1[1]:0:-12}
  • macではpfx=${r1[$i]:0:$((${#r1[$i]}))-12}

# シングル
idx="~/db/index/hisat2_idx" # リファレンスインデックス
r1=(*_R1.fastq.gz) # R1ファイル名の配列
nsfx=12 # 拡張子文字数
for i in ${r1[@]}; do
 pfx=${i:0:$((${#i}))-$nsfx}
 hisat2 -x $idx -U ${i} -p 4 2>> log.txt | samtools sort -@ 4 > ${pfx}.sort.bam
done

# ペアエンド
idx="~/db/index/hisat2_idx" # リファレンスインデックス
r1=(*_R1.fastq.gz) # R1ファイル名の配列
r2=(*_R2.fastq.gz) # R2ファイル名の配列
n=${#r1[@]} # ファイル数

for i in `seq 0 $((--n))`; do
 pfx=${r1[$i]:0:$((${#r1[$i]}))-12}
 hisat2 -x $idx -1 ${r1[$i]} -2 ${r2[$i]} -p 4 2>> log.txt | samtools sort -@ 4 > ${pfx}.sort.bam
done

##  [1] "100000 reads; of these:"                                                
##  [2] "  100000 (100.00%) were paired; of these:"                              
##  [3] "    27240 (27.24%) aligned concordantly 0 times"                        
##  [4] "    66390 (66.39%) aligned concordantly exactly 1 time"                 
##  [5] "    6370 (6.37%) aligned concordantly >1 times"                         
##  [6] "    ----"                                                               
##  [7] "    27240 pairs aligned concordantly 0 times; of these:"                
##  [8] "      1318 (4.84%) aligned discordantly 1 time"                         
##  [9] "    ----"                                                               
## [10] "    25922 pairs aligned 0 times concordantly or discordantly; of these:"
## [11] "      51844 mates make up the pairs; of these:"                         
## [12] "        29890 (57.65%) aligned 0 times"                                 
## [13] "        16949 (32.69%) aligned exactly 1 time"                          
## [14] "        5005 (9.65%) aligned >1 times"                                  
## [15] "85.06% overall alignment rate"

2 samtools

---

sam形式のアラインメントファイルを扱うためのユーティリティ

• bamに変換samtools view
  • b : 出力がbam
  • S : 入力がsamの場合(旧バージョンのみ,現在は自動判定)
• リードをポジション順にソート. samtools sort
• bamインデックス(*.bai)を作成samtools index
• マップ率の値: samtools flagstat
• samtools idxstat
• デプスsamtools depth
• bamファイルをマージ samtools merge
• PCR duplicateを除く samtools rmdup

# bam -> sam

# sam -> bam  
samtools view -b seq.sam > seq.bam

# bamをsort
samtools sort seq.bam > seq.sort.bam

# sortコマンドに直接samファイルを指定
samtools view -b seq.sam -@ 4 | samtools sort > seq.sort.bam
samtools sort -@ 4 seq.sam > seq.sort.bam

# unmapped / mapped bam
samtools view -b -f 4 file.bam > unmapped.bam
samtools view -b -F 4 file.bam > mapped.bam

# bam index
samtools index seq.sort.bam -@ 4

# mapping, bam変換, sortを連続して実行
bowtie2 -x ref seq.fastq -p4 | samtools view -b - | samtools sort > seq.sort.bam
bowtie2 -x ref seq.fastq -p4 | samtools sort -@ 4 > seq.sort.bam

# samtools flagstat
samtools flagstat seq.sort.bam > seq.flgst

# samtools merge
samtools merge merged.bam 1.bam 2.bam

2.1 ソート済bamか否かを確認

samtools view -H seq.bam | grep "SO:"
samtools view -H seq.sorted.bam | grep "SO:"

## @HD  VN:1.0  SO:unsorted
## @HD  VN:1.0  SO:coordinate

2.2 マップ率 samtools flagstat, bamtools

• properly paired: 同じクロモソームに適切な向きと位置でアライメントされたペアリードの数

samtools flagstat seq.bam 

##  [1] "213600 + 0 in total (QC-passed reads + QC-failed reads)"
##  [2] "13600 + 0 secondary"                                    
##  [3] "0 + 0 supplementary"                                    
##  [4] "0 + 0 duplicates"                                       
##  [5] "183710 + 0 mapped (86.01% : N/A)"                       
##  [6] "200000 + 0 paired in sequencing"                        
##  [7] "100000 + 0 read1"                                       
##  [8] "100000 + 0 read2"                                       
##  [9] "145520 + 0 properly paired (72.76% : N/A)"              
## [10] "150132 + 0 with itself and mate mapped"                 
## [11] "19978 + 0 singletons (9.99% : N/A)"                     
## [12] "3966 + 0 with mate mapped to a different chr"           
## [13] "3133 + 0 with mate mapped to a different chr (mapQ>=5)"

2.3 リファレンス各位置にマップされたリード数: samtools depth

samtools depth -a seq.sort.bam | gzip > seq.dep.gz

2.4 bamからランダムにリードを取り出す

2.4.1 bamtools random

2.4.2 sort -R | head -n

samtools view input.bam | sort -R | head -n 100 | samtools view -b > output.bam

## [E::hts_open_format] Failed to open file input.bam
## samtools view: failed to open "input.bam" for reading: No such file or directory

2.4.3 shuf -n

shuf -n N input > output

3 発現解析

---

3.1 stringtie

RNA-seqの場合アラインメントデータ(bam)を元にして, transcriptを構築する(gtfが作られる). 構築されたtranscriptごとに定量することがができる.

• primary GTFの作成
  • -G: ガイドgff
  • -m: 最小トランスクリプト長 (default: 200)
• gtfファイルのマージ(gtfファイルのリストを入力とする.)
• マージされたgtfファイルをガイドとして定量する
  • -e: リードカウントを-Gで指定されたgtfの領域に制限.
  • -A: カウントデータの出力ファイル
  • -b:カバレッジデータを含むBallgown入力テーブルファイル（* .ctab）を出力

# bamファイルからgtfの作成(サンプルごとに作成)
guide=~/db/index/hisat2_idx/ref.gff # ガイドgff
sfx=.sort.bam # 拡張子
nsfx=${#sfx} # 拡張子文字数
for f in *.sort.bam; do 
  pfx=${f:0:$((${#f}))-$nsfx}
  stringtie -p 4 ${f} -G ${guide} -o ${pfx}.gtf
done

# gtfファイル(フルパス)のリストを作成
greadlink -f `ls *.gtf` > gtf.list.txt # mac
readlink -f `ls *.gtf` > gtf.list.txt # ubuntu 

# gtfのマージ gtfファイルパスを記述したテキストファイル(gff.list.txt)を入力
stringtie --merge -G guide_gff -o merged.gtf gtf.list.txt

# マージされたgtfファイルをガイド(-G)として指定して, その領域に絞って定量する
for f in *.sort.bam; do 
  pfx=${f:0:$((${#f}))-$nsfx}
  stringtie -p 4 ${f} -e -A ${pfx}.cnt.tab -b res_ballgown -G merged.gtf -o ${pfx}.gtf
done

3.2 rsem　

リードアラインメントと定量化を一括して実行できる. alignerをbowtie2, STARから選択する. transcriptとgeneのリストファイルを用いてリファレンスインデックスを作成すれば, それぞれ別個にカウントすることができる.

• rsem-calculate-expression
  • 標準出力にコマンドやログが, 標準エラー出力にマップ率の結果が出力される.
  • --strand-specific : bowtie2の -- norc オプションがつく
  • -p:スレッド数指定
  • --bowtie2
  • --no-bam-output bamファイルを作らない.
  • --sort-bam-by-coordinate ソートされたbamを作る.
  • --strandednes <none|forward|reverse> norc’ / ’

# index作成
rsem-prepare-reference --bowtie2 ref.fasta ref

# transcriptとgeneのリストファイルを用意して, index作成
rsem-prepare-reference --bowtie2 --transcript-to-gene-map list.txt ref.fasta ref

# single
rsem-calculate-expression --bowtie2 -p 4 seq_R1.fastq.gz ref out
rsem-calculate-expression --bowtie2 -p 4 --sort-bam-by-coordinate seq_R1.fastq.gz ref out
rsem-calculate-expression --bowtie2 -p 4 seq_R1.fastq.gz ref out >log 2> err.txt

# pair
rsem-calculate-expression --bowtie2 --sort-bam-by-coordinate --paired-end seq_R1.fastq.gz seq_R2.fastq.gz ref out_dir

3.2.1 リストファイルの作成

• 遺伝子ID トランスクリプトIDのようなファイルを作成する

# 例
x <- "XM_003495724.4 PREDICTED: Cricetulus griseus heat shock protein family B (small) member 11 (Hspb11), transcript variant X1, mRNA"
regmatches(x, gregexpr("\\(.*?\\)", x))

## [[1]]
## [1] "(small)"  "(Hspb11)"

ref <- Biostrings::readDNAStringSet("~/pub/dat/191007_sumitomo/dat/shina/GCF_000223135.1_CriGri_1.0_rna_chrMT_AB.fa")

add <- names(ref)[46718:46721]
href <- names(ref)[-46718:-46721]

gid <- unlist(lapply(href, function(x){
  v <- unlist(regmatches(x, gregexpr("\\(.*?\\)", x)))
  sub("\\)", "", sub("\\(", "", v[length(v)]))
  }))
tid <- sapply(strsplit(href, " "), "[", 1)

tmp <- data.frame(x = gid, y = tid, stringsAsFactors = F)

3.3 salmon

mapping based modeとalignment-based modeがある. alignment-based modeでは, 好みのmapperで出力したbamを入力として使える.
- salmon index -t ref.fa -i ref - インデックスディレクトリができる. - インデックスディレクトリ中のduplicate_clusters.tsvに重複したトランスクリプトのリスト

• salmon quant --help-reads, salmon quant --help-alignment
• -a bamファイル
• -o で出力ディレクトリ指定.
  • quant.sfカウントデータが中に作られる.
• -l ライブラリタイプ. library typeについてFragment Library Types
  • A : bamインプットの場合自動判定
  • SR: ストランデッド. read1がリバースストランド
  • ISR: ペアエンド内向き. ストランデッド. read1がリバースストランド.
• --validateMappings

# index作成 -kでkmer指定[default:31]
salmon index -p 4 -t ref.fa -i ref
salmon index -p 4 -k 31 -t ref.fa -i ref

# single
salmon quant -p 4 -i ref -l A -r seq_R1.fastq.gz --validateMappings -o res

# pair Quantifying in mapping-based mode
salmon quant -p 4 -i ref -l A -1 seq_R1.fastq.gz -2 seq_R2.fastq.gz - o res

idx=~/db/index/salmon_idx/mfq_tr3_cdh # リファレンスインデックス
r1=(*.fastq.gz) # R1ファイル名の配列
nsfx=9 # 拡張子文字数
mkdir ../output # 出力ディレクトリ

for i in ${r1[@]}; do
 pfx=${i:0:$((${#i}))-$nsfx}
 # echo $i; echo $pfx
 salmon quant -p 4 -i ${idx} -l A -r ${i} --validateMappings -o ../output/${pfx}
done

3.4 kallisto

• --single Quantify single-end reads
• -l 64.3 平均フラグメント長[default: インプットから推定される.]
• -s 2.2
• -t thread数(bootstrapを実行する場合)

kallisto index -i ref ref.fa
kallisto quant -i ref -o <output_dir> -b <bootstrap_num> --single -l <fragment_length>

4 補足: ベースコール (*.bcl)ファイルからfastqファイルを生成する.
以下は, TruSeq CD Indexesを使用する場合.

• indexの配列を確認するillumina-adapter-sequences-1000000002694-07.pdf.
  • Illumina Adapter Sequences (1000000002694 v10)を探してindex配列を確認
  • TruSeqの場合はTruSeq DNA and RNA CD(Combinatorial dual) Indexes を見る.
• indexのレイアウトを確認する.index-adapters-pooling-guide-1000000041074-06
  • Index Adapters Pooling Guideを探す.
  • DAP or RAP Dual-Indexed Layoutを見ること

4.1 SampleSheet.csvファイルを作成

• エクセルで編集すると^MがついてくるのでRで編集して書き出す(macで作業している場合).
• write.csvではなくwrite.table(sep=“,”, col.names=F)にする

4.2 bcl2fastq

• runディレクトリ内にSampleSheet.csvをあらかじめ置いておく.

# bcl2fastq実行
bcl2fastq --run-folder-dir run_dir -p 12 --output-dir run_dir/fastq_files --no-lane-splitting

# コマンド省略
bcl2fastq -R run_dir -p 12 -o run_dir/fastq_files --no-lane-splitting

# runディレクトリ内でbcl2fastqを実行
cd ./run_dir
bcl2fastq -p 12 -o ./fastq_files --no-lane-splitting

# fastqファイルをサンプルディレクトリから移動する
mkdir fastq
smp=(`ls`)
for d in ${smp[@]}; do
mv ./$d/*.fastq.gz ./fastq
done

5 環境

sessionInfo()

## R version 3.5.2 (2018-12-20)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] ja_JP.UTF-8/ja_JP.UTF-8/ja_JP.UTF-8/C/ja_JP.UTF-8/ja_JP.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] dplyr_0.8.3
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.3          XVector_0.22.0      knitr_1.24         
##  [4] magrittr_1.5        zlibbioc_1.28.0     IRanges_2.16.0     
##  [7] BiocGenerics_0.28.0 hms_0.5.1           tidyselect_0.2.5   
## [10] R6_2.4.1            rlang_0.4.2         stringr_1.4.0      
## [13] tools_3.5.2         parallel_3.5.2      xfun_0.9           
## [16] htmltools_0.3.6     yaml_2.2.0          digest_0.6.23      
## [19] assertthat_0.2.1    tibble_2.1.3        crayon_1.3.4       
## [22] purrr_0.3.3         readr_1.3.1         S4Vectors_0.20.1   
## [25] vctrs_0.2.1         zeallot_0.1.0       glue_1.3.1         
## [28] evaluate_0.14       rmarkdown_1.15      stringi_1.4.3      
## [31] compiler_3.5.2      pillar_1.4.2        Biostrings_2.50.2  
## [34] backports_1.1.5     stats4_3.5.2        pkgconfig_2.0.3