3.2. STAR aligner로 mapping 하기

우선 참조유전체(reference genome)이 있어야 한다.

그리고 paired end ILLUMINA sequencing fastq를 사용하는 이유는 요새는 이게 가장 싸다.
(왜냐하면 NovaSeq으로 돌려야 bp당 sequencing 비용이 가장 저렴한데, 한번 flowcell을 돌릴때 같은 length의 single-end or paired-end data를 generation 해야한다. NovaSeq은 WGS 기준 30X로 30명 분을 생산한다고 하기 때문에, 다 같은 형태의 data - 같은 single or pair end, 그리고 150 bp or 100 bp- 형태로 같게 해야 한다.)

그렇기 때문에 나의 안내도 마찬가지로 paired-end mapping으로 기록하였다. 

준비물: reference fasta (indexed)

STAR --runThreadN [thread] \
    --genomeDir [reference genome directory]\
    --outSAMtype BAM Unsorted \
    --quantMode GeneCounts \
    --outFileNamePrefix [prefix] \
    --readFilesCommand zcat \
    --readFilesIn [input fastq pair 1] [input fastq pair 2]

--runThreadN에서 thread: 계산의 병렬 정도를 조절 가능함 (일반적으로 CPU core 하나는 hyper thread에 따라 2 threads이기 떄문에 가능한 계산 자원을 파악하라)
--outSAMtype에서 BAM unsorted가 속도 측면에서 유리해서 이용하였다.
--quantMode에서 GeneCounts 옵션을 주면 HTseq으로 gene counts 한 결과와 같은 결과를 얻을 수 있다.
--readFilesCommand에서 zcat을 쓰는 이유는 fastq.gz 를 input으로 사용하기 때문이다.

 

이렇게 수행하게 되면 HTseq의 결과와 같은 gene counts를 얻을 수 있고, genome viewer에서 read를 볼 수 있는 BAM을 확보 할 수 있다. 

물론 genome viewer에서 보려면 다음과 같이 sorting 과정과 indexing 과정이 필요하다. 

samtools sort -m 10G -@ [thread] -o out [sorted bam] [input bam]
Samtools index [sorted bam]

-m: 10G 이유는 일반적으로 1개의 thread가 10G를 쓴다고 해서 준 수치이다. memory가 클수록 더 효율적으로 빠르게 분석이 된다.
단지 여기서 thread * m 으로 메모리를 먹기 때문에 이를 고려해야 한다.