Fastq to consensus sequences

 

 

Bacteria genome sequencing을 Pacbio로 했다면 필요없는 과정이겠지만, Miseq이나 Hiseq으로 데이터를 만들었고 이걸 하나의 sequence로 합치고자 한다면 아래와 같은 과정이 필요합니다. Denovo sequencing을 하려는 경우는 해당되지 않아요. Reference 에 mapping하고 하나의 sequence를 만드는 과정입니다. Sample이 하나의 chromosome으로 되어있지않고 Plasmid가 하나 또는 여러개가 더 있을 경우 아래 있는 방법으로 하면 plasmid sequence는 놓칠 수 있어요.

 

- SAMtools와 BCFtools는 버젼을 같은 것으로 하는게 좋습니다. 둘 다 한꺼번에 새버젼으로 새로 설치를 하는게 좋습니다. 버젼이 다르면 오류가 생길 수 있습니다. 

- 분석 때 생기는 오류들을 다음 글을 참고해서 수정하면서 만들었습니다.

 

 

1. https://solution4u.tistory.com/10

3.1. BWA를 통한 FASTQ의 맵핑 (BWA) 진행하기

2. http://blog.genoglobe.com/2015/09/next-generation-sequencing-variant.html

[Next-Generation Sequencing] Variant calling 작업의 오류 해결과 Bowtie 2

 

Summary

명령어 summary 입니다. 자세한 해석이 필요없고, 원래 해봤는데 순서나 명령어 일부가 기억이 안나신 분들은 여기를 보시면 됩니다. 

bwa index –a bwtsw reference.fasta
bwa mem –t 10 ref.fasta sample_1.r1.fastq sample_1.r2.fastq > aln-pe.sam
samtools view –bS aln-pe.sam > aln-pe.bam
samtools sort aln-pe.bam > aln-pe.sorted.bam
samtools index aln-pe.sorted.bam
bcftools mpileup –Ou -d 8000 –f reference.fasta aln-pe.sorted.bam | bcftools call –mv –Oz –o calls.vcf.gz
bcftools index calls.vcf.gz
cat reference.fasta | bcftools consensus calls.vcf.gz > consensus.fasta

---

 

1. Reference genome에 mapping - 사용 프로그램(버전) = BWA (0.7.12-r1039)

 

1.1 Reference sequence index file 만들기

bwa index [-p prefix] [-a algotype] xx.fasta  # 예시

• a 옵션은 두 가지(is, bwtsw) 있습니다. is는 database가 2GB 이상일 때 작동 안 한다고 합니다. Database가 2GB라는게 무슨말일까요....일단 bwtsw를 써서 진행했습니다.

• -p 옵션은 레퍼런스 이름을 적는 것 같습니다(same as db filename). – 이거가 파일 이름이랑 같아야 하는 듯. 이걸 –p EDL933으로 하면 index 관련 파일들의 이름이 EDL933.sa 등으로 만들어짐. 이때 reference 파일 이름이 EDL933.fasta 처럼 되어 있어야 함. -p는 일단 사용할 필요가 없는 것 같습니다. 언제 필요하긴 할 것 같은데 없어도 일단 index 파일이 잘 만들어집니다.

 

그래서 실제로 사용할 명령어는:

bwa index –a bwtsw reference.fasta

 

1.2. Reference sequence에 mapping 하기

bwa mem –M –R @RG\tID:HWI\tSM:4\tPL:ILLUMINA\tLB:miseq –t 10 EDL933.fna file1.fastq.qz file2.fastq.qz

이 명령어는 여기 보고 만들었는데요...안됩니다...뭔가 오타가 있는 것 같은데..다른 글을 참고해서 다음 과정을 진행했습니다. 

bwa mem –t 10 ref.fasta sample_1.r1.fastq sample_1.r2.fastq > aln-pe.sam

–t 옵션은 mapping 에 쓸 cpu 개수를 정해주면 됩니다.

pe는 paired end (하나의 sample을 paire end 방식으로 한 것)로 염기서열을 얻어서 이렇게 적었습니다. 이거는 이렇게 안 적어도 무방하나 나중에 만들어진 sam 파일을 single end 파일로 만든건지 paired end 파일로 만든건지 알기 쉽게 분류하기 위해 이렇게 적었습니다.

 

2. BAM file 만들기

samtools view –bS aln-pe.sam > aln-pe.bam

 

3. Sorted bam file 만들기 (중요!)

samtools sort aln-pe.bam > aln-pe.sorted.bam

bam file index만들 때 bam 파일이 필요한데요, 그 전에 bam 파일을 sorting 한 sorted bam 파일을 만들어야 consensus sequence 만들때 오류 없이 잘 진행됩니다.

 

4. Bam file index 만들기 (bam.bai)

samtools index aln-pe.sorted.bam

이렇게 하면 bam.bai 파일이 만들어집니다.

 

5. Reference file index 만들기 (fasta.fai)

  - "1" 에 설명되어 있습니다. reference sequence의 인덱스(.fai) 파일이 없는 경우 진행하시면 되고 이미 있다면 넘어가면 됩니다.

 

6. Consensus sequence 만들기

 

6.1 VCF 파일 만들기

bcftools mpileup –Ou -d 8000 –f reference.fasta aln-pe.sorted.bam | bcftools call –mv –Oz –o calls.vcf.gz

- 원래 –d의 default 값은 250 입니다. File당 read 수라는데...세어보지 않았지만 대충 read수는 그 파일에 250일거보다 많은 것 같았고, 구글링하니 사람들이 bacteria 다룰 때 –d 값을 8000으로 하는 것 같아서 이렇게 했습니다.

 

6.2. VCF 파일 index 하기

bcftools index calls.vcf.gz

 

6.3. consensus fasta file 만들기

cat reference.fasta | bcftools consensus calls.vcf.gz > consensus.fasta

 

이렇게 하면 하나로 연결된 DNA 염기서열 데이터가 만들어집니다.

 

- 진행하다보면 오류가 생기실텐데요, 이럴 경우 명령어를 제대로 쳤는지 봐야합니다. 6.1에 있는 명령어 중 – mv의 경우 - 이상하게도 꼭 -my로 타이핑 되어있습니다. 똑같은걸 몇 번이나 반복해서 실수 한 적이 있네요...

- 한번 더 말씀드리지만 bcftools, samtools http://www.htslib.org/download/등은 업데이트 버전 받을 때 두 개를 같이 받아야 합니다. 호환되는 정보, 명령어들이 조금씩 달라지는 것 같습니다...

SAMtools/BCFtools/HTSlib - Downloads

여기 보면 HTSlib 가 있는데요 이거는 안 받아도 consensus sequence가 만들어집니다. 이건 뭐하는건지 아직 모르겠네요...