====== 일루미나 데이터의 QC와 기본 전처리 ======
Genome sequencing의 기본 목표는 (1) de novo assembly를 통한 새로운 유전체 서열의 구성, 또는 (2) re-sequencing을 통한 변이의 예측이다. Reference genome sequence가 필요한 것은 후자의 경우이다. 무엇을 목표로 하든 일루미나 장비에서 생산된 시퀀싱 raw data를 평가(데이터 자체를 변경하지 않는 탐색적 분석)하고 필요하다면 trimming 등의 적극적인 조작을 해야 한다.

===== Sequence QC =====

가장 대중적인 QC 프로그램은 [[https://www.bioinformatics.babraham.ac.uk/projects/fastqc/|FastQC]]이다. 이 프로그램은 read data에 대하여 아무런 조작을 가하지 않고 QC 관련 plot만 생성한다. 명령행 혹은 GUI 모드로 전부 실행 가능하다. 단순히 read의 수, 평균 길이, 총 염기수 등을 집계하려면 [[https://github.com/dib-lab/khmer|khmer]] 패키지의 readstats.py를 사용하면 된다. 어댑터 제거 및 quality trimming을 실시하려면 [[https://github.com/timflutre/trimmomatic|trimmomatic]]을 사용한다.

===== 일반적인 서열 데이터 조작 =====

[[https://emboss.sourceforge.net/|EMBOSS package]]에는 다수의 기본적인 생명정보학 분석 프로그램이 포함되어 있다. wossname <keyword>를 입력하면 short description에 keyword가 포함된 개별 프로그램의 목록이 나온다. tfm <program>을 입력하면 각 프로그램에 대한 상세한 매뉴얼이 출력된다. 서열 파일의 입력과 출력 및 포맷 변환 등을 수행하는 seqret에 대해서 알아보자.

  # Multiple fasta file을 개별 서열 파일로 분리
  # -auto 옵션은 입력하지 않은 값에 대한 프롬프트가 나오지 않게 함
  $ seqret -ossinlge -auto contigs.fasta
  # 여러 서열로 구성된 GenBank를 개별 파일로 분리하되 feature를 전달
  $ seqret –feature –ossingle –os format genbank –auto INPUT.gbk
  # 주어진 단일서열 파일의 5 염기와 끝에서 5번째 염기까지를 추출하여 reverse complementary 형태로 전환한 뒤 별도의 파일로 저장
  $ seqret –sbegin 5 –send -5 –sreverse INFILE.fasta OUTFILE.fasta
  # 다중서열 파일에서 특정 ID의 서열을 별도의 파일로 추출
  $ seqret all.fasta:seq_id -auto

[[https://jgi.doe.gov/data-and-tools/software-tools/bbtools/|BBTool]]은 DNA 및 RNA 서열을 위한 다양하고 방대한 분석 도구의 모음이다. 예를 들어 reformat.sh는 read length 분포의 히스토그램을 계산하고, stats.sh는 assembly statistics를 산출한다. BBTool 패키지에 포함된 bbnorm.sh는 k-mer 기반의 normalization을 수행하는 용도로 만들어진 것인데, 오류 교정이나 k-mer frequency plot을 그리는 데에도 쓸 수 있다. 자세한 설명은 SEQanswers의 [[https://www.seqanswers.com/forum/bioinformatics/bioinformatics-aa/44377-introducing-bbnorm-a-read-normalization-and-error-correction-tool?t=49763|Introducing BBNorm, a read normalization and error-correction tool]]을 참조하기 바란다.

[[https://github.com/lh3/seqtk|seqtk]]는 FAST/Q format file의 조작을 위한 도구이다. 일례로 다음은 대용량의 fastq 파일에서 1만개의 서열을 subsampling하는 방법이다. Pair 상태를 유지하려면 같은 random seed 값(-s100)을 사용해야 한다.

  $ seqtk sample –s100 read_1.fq 10000 > sub_1.fq
  $ seqtk sample –s100 read_2.fq 10000 > sub_2.fq

일루미나 데이터를 위한 QC 및 전처리 도구로서 [[https://github.com/jts/sga/wiki/Preqc|SGA preqc]]라는 것도 있다. Preqc는 de novo assembler의 일종인 [[https://github.com/jts/sga|SGA(String Graph Assembler)]]의 부속 스크립트인데, sequence coverage, per-based error rate, genome size, heterozygosity 및 repeat content 등 다양한 지표를 추정하는 역할을 한다. 실행시간이 다소 오래 걸린다는 것이 단점이다. 

  $ sga preprocess --pe-mode 1 reads_R1.fastq reads_R2.fastq > mygenome.fastq
  $ sga index -a ropebwt --no-reverse -t 8 mygenome.fastq
  $ sga preqc -t 8 mygenome.fastq > mygenome.preqc
  $ sga-preqc-report.py mygenome.preqc sga/src/examples/*.preqc
  
[[https://github.com/TGAC/KAT|KAT]]는 k-mer 분석을 통해서 NGS 데이터셋과 유전체 조립물에 대한 QC를 실시하는 도구이다. 그래픽을 수록한 보고자료 생성, 여러 데이터셋에 대한 k-mer count 비교, 필터링 등이 한번에 실시되므로 매우 편리하다. K-mer 길이를 –m <num>으로 지정하지 않으면 27을 기본 수치로 사용한다. Kat hist는 jellyfish와 유사한 k-mer 분포 히스토그램을, kat gcp는 각 k-mer에 대하여 GC count를 계산한 매트릭스를 작성한다. 두 명령어 모두 png 파일 형태의 diagnostic plot을 출력한다.

  # OUTPUT_PREFIX를 지정하지 않으면 'kat.command'가 사용됨
  $ kat hist –t 8 –o OUTPUT_PREFIX reads_R1.fastq reads_R2.fastq
  $ kat gcp -t 8 reads_R1.fastq reads_R2.fastq
  $ kat comp 'reads_R1.fastq reads_R2.fastq' assembly.fa
  
[[https://expressionanalysis.github.io/ea-utils/|ea-utils]]는 fastq 파일로부터 barcode demultiplexing, paired end joining, adapter trimming 등을 수행하는 유틸리티 모음이다. 앞에서 설명한 여러 종류의 분석 또는 조작을 수행하는 프로그램들이 한 패키지 안에 포함되어 있으므로 적극 활용을 권장한다.

  $ fastq-stats FILE.fastq # 표준 출력으로 sequence statistics 수치를 출력함

===== SSU rRNA 서열의 재구성 ======
[[https://github.com/HRGV/phyloFlash|phyloFlash]]는 일루미나 데이터셋으로부터 SSU rRNA 서열을 재구성하여 동정하는 도구이다. Silva reference database에 대한 read mapping을 하여 해당되는 read를 회수한 뒤 targeted assembly를 거쳐 전장 SSU rRNA 서열을 만들어낸다. 원래는 메타게놈 혹은 메타트랜스크립톰 read의 분석을 위해 만들어진 것이지만 단일 유전체의 시퀀싱 read에 이를 적용하면 시퀀싱 라이브러리의 오염 여부를 점검할 수도 있다. 결과물은 단순 텍스트 및 HTML 형식의 리포트로 제공된다. 아래에 보인 실행 사례에서 모든 결과 파일은 LIB이라는 접두사를 갖는다.

  $ phyloFlash.pl -lib LIB -read1 reads_F.fq.gz -read2 reads_R.fq.gz -CPUs 16
  $ phyloFlash.pl -lib LIB -CPUs 16 -read1 reads_FR.fq.gz # interleaved file