Differences

This shows you the differences between two versions of the page.

--- bioinfo:whole-genome_alignment [2023/06/28 08:48] – [minimap2의 활용] hyjeong
+++ bioinfo:whole-genome_alignment [2023/06/29 10:05] (current) – [dnadiff 결과로부터 균주 특이적 영역 찾기] hyjeong
@@ Line 33: / Line 33: @@
   * out.unqry   - Unaligned query sequence IDs and lengths
-리포트 파일(out.report)의 내용은 다음과 같다. Refenece와 query 전부 complete genome sequence이며, reference는 16 kb의 플라스미드가 있다. 리포트 파일 내용의 설명은 [[https://github.com/marbl/MUMmer3/blob/master/docs/dnadiff.README|dnadiff REAME]] 파일을 참조하라.
+리포트 파일(out.report)의 내용은 다음과 같다. Refenece(GCF_000063585.1)와 query(GCF_000017045.1) 전부 complete genome sequence이며, reference는 16 kb의 플라스미드가 있다. 리포트 파일 내용의 설명은 [[https://github.com/marbl/MUMmer3/blob/master/docs/dnadiff.README|dnadiff REAME]] 파일을 참조하라.
   /home/hyjeong/temp_20200624/ATCC_3502.fna /home/hyjeong/temp_20200624/Hall.fna
@@ Line 52: / Line 52: @@
   UnalignedSeqs               0(0.00%)             0(0.00%)
   <이하 생략>
+SNP 파일을 VCF로 전환하려면 [[https://github.com/MatteoSchiavinato/all2vcf|all2vcf]] 유틸리티의 mummer 명령을 이용하라.
 두 염기서열 FASTA file의 identity만을 계산하고 싶다면 identity.sh 스트립트를 사용하라. 내가 만든 스크립트였던가, 혹은 어디서 퍼 온 것인가? 잘 기억이 나지 않는다.
@@ Line 143: / Line 145: @@
   $ seqtk seq -L 150 query-specific.fa > filtered.fa
-유전체 염기서열의 특정 구간 정보를 BED 파일로 확보하였다면 여기에 어떤 기능을 갖는 유전자가 위치하고 있는지를 알고 싶을 것이다. GenBank 파일에서 유전자 정보를 수록한 BED 파일을 만든 뒤, 이를 위에서 만든 query-specific.bed(.filt)와 함께 bedtools intersect 명령어에 인수로 제공하면 된다. 첫 번째 BED 파일에는 유전자의 locus tag과 기능 정보를 4번째 컬럼부터 삽입해 두면 최종 결과 파일에도 이 정보가 상속되므로 후속 작업에 매우 편리하다. /data/BinScript/gbkInfo.pl 결과물을 활용하면 첫 번째 BED 파일을 만들 수 있다. 두 BED 파일의 경계에 걸치는 유전자를 처리하기 위하여 [[https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html|bedtools intersect]] 명령 실행 시 '-f 0.3' 옵션을 주었다. 이는 유전자 정보를 담고 있는 첫 번째 BED 파일의 정보에 대하여 최소 30%가 겹쳐야 함을 의미한다. 이는 단지 하나의 사례이므로 적절히 바꾸어도 된다. 예를 들어 '-f 1.0'으로 설정하면 specific region 안에 전 영역이 완전히 포함되는 유전자들만 추출될 것이다.
+유전체 염기서열의 특정 구간 정보를 BED 파일로 확보하였다면 여기에 어떤 기능을 갖는 유전자가 위치하고 있는지를 알고 싶을 것이다. GenBank 파일에서 유전자 정보를 수록한 BED 파일을 만든 뒤, 이를 위에서 만든 query-specific.bed(.filt)와 함께 bedtools intersect 명령어에 인수로 제공하면 된다. 첫 번째 BED 파일에는 유전자의 locus tag과 기능 정보를 4번째 컬럼부터 삽입해 두면 최종 결과 파일에도 이 정보가 상속되므로 후속 작업에 매우 편리하다. [[:custom_scripts#gbkinfopl|gbkInfo.pl]]을 활용하면 첫 번째 BED 파일을 만들 수 있다. 두 BED 파일의 경계에 걸치는 유전자를 처리하기 위하여 [[https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html|bedtools intersect]] 명령 실행 시 '-f 0.3' 옵션을 주었다. 이는 유전자 정보를 담고 있는 첫 번째 BED 파일의 정보에 대하여 최소 30%가 겹쳐야 함을 의미한다. 이는 단지 하나의 사례이므로 적절히 바꾸어도 된다. 예를 들어 '-f 1.0'으로 설정하면 specific region 안에 전 영역이 완전히 포함되는 유전자들만 추출될 것이다.
   $ gbkInfo.pl Hall_A.gbk # Hall_A.gbk.txt 파일 생성(query genome)
@@ Line 161: / Line 163: @@
 ===== minimap2의 활용 =====
-minimap([[https://academic.oup.com/bioinformatics/article/32/14/2103/1742895|논문]])은 원래 long read와 reference의 빠르고도 대략적인 alignment를 위해서 만들어진 프로그램이다(심지어 20%의 차이까지 수용하는!). 현재 쓰이는 [[https://github.com/lh3/minimap2|minimap2]]는 genome과 spliced nucleotide의 정렬 목적으로 쓰일 수 있게 더욱 발전하였다. minimap2의 여러 사용례 중에서 미생물 유전체를 1:1로 비교하여 SNP를 추출하는 방법을 정리해 보았다. 이를 더욱 확장하면 정렬이 이루어지지 않는 영역 정보를 추출하는 것도 가능할 것이다.
+이 항목은 [[bioinfo:참조_서열에_대한_매핑_reference_mapping_및_시각화|참조 서열에 대한 매핑(reference mapping) 및 시각화]]에서 다루는 것이 원칙적으로 맞을 것이다. 그러나 reference에 대하여 다수의 짧은 read를 붙이는 것이 아니라 이미 조립이 완료된 긴 염기서열을 1:1로 비교하는 것을 주된 내용으로 하므로, whole-genome alignment 항목에서 설명하는 것이 더 낫다고 생각한다.
-[[https://github.com/lh3/minimap2/tree/master/misc|minimap misc]] 문서의 [[https://github.com/lh3/minimap2/tree/master/misc#calling-variants-from-haploid-assemblies|Calling
+minimap([[https://academic.oup.com/bioinformatics/article/32/14/2103/1742895|논문]])은 원래 long read와 reference의 빠르고도 대략적인 alignment를 위해서 만들어진 프로그램이다(심지어 20%의 차이까지 수용하는!). 이는 BWA의 개발자인 [[http://liheng.org/|Heng Li]]가 개발하였다. 현재 쓰이는 [[https://github.com/lh3/minimap2|minimap2]]([[https://academic.oup.com/bioinformatics/article/34/18/3094/4994778|2018년도 논문]])는 genome과 spliced nucleotide의 정렬 목적으로 쓰일 수 있게 더욱 발전하였다. Heng Li의 블로그에서 minimap2로의 전환 및 BWA의 미래에 관하여 쓴 글([[https://lh3.github.io/2018/04/02/minimap2-and-the-future-of-bwa|링크]])을 참조하라. 또한 공식 문서의 [[https://github.com/lh3/minimap2#getting-started|Getting Started]]를 통해서 minimap2가 여러 상황에서 어떻게 쓰이는지 감상해 보기 바란다.
-variants from haploid assemblies]] 항목을 주로 참조하였다.
+년도 논문에서는 minimap2의 특성을 다음과 같이 소개하였다.
+<color #00a2e8>Minimap2 is a versatile mapper and pairwise aligner for nucleotide sequences. It works with short reads, assembly contigs and long noisy genomic and RNA-seq reads, and can be used as a read mapper, long-read overlapper or a full-genome aligner. Minimap2 is also accurate and efficient, often outperforming other domain-specific alignment tools in terms of both speed and accuracy.</color>
+minimap2의 여러 사용례 중에서 미생물 유전체를 1:1로 비교하여 SNP를 추출하는 방법을 정리해 보았다. 이를 더욱 확장하면 정렬이 이루어지지 않는 영역 정보를 추출하는 것도 가능할 것이다.
+[[https://github.com/lh3/minimap2/tree/master/misc|minimap2 misc]] 문서의 [[https://github.com/lh3/minimap2/tree/master/misc#calling-variants-from-haploid-assemblies|Calling
+variants from haploid assemblies]] 항목을 주로 참조하였다. minimap2의 여러 옵션에 대해서는 따로 공부를 하는 것이 좋을 것 같다. PAF(a Pairwise mApping Format) 파일의 설명은 [[https://github.com/lh3/miniasm/blob/master/PAF.md|여기]]를 참조하라. paftools.js는 minimap2과 같이 설치되었으므로 어디에 있는지 잘 찾아보기 바란다. 아래의 실행 사례는 위에서 dnadiff 분석의 샘플로 사용했던 보툴리늄 균주 2종의 유전체를 사용한 것이다.
   # keeping this file is recommended; --cs required!
@@ Line 171: / Line 181: @@
   $ sort -k6,6 -k8,8n asm.paf > asm.srt.paf
   $ k8 paftools.js call asm.srt.paf > asm.var.txt
+  3686527 reference bases covered by exactly one contig
+substitutions; ts/tv = 3.671
+1bp deletions
+1bp insertions
+2bp deletions
+2bp insertions
+[3,50) deletions
+[3,50) insertions
+[50,1000) deletions
+[50,1000) insertions
+>=1000 deletions
+>=1000 insertions
+  $ head -n 5 asm.var.txt
+  R       NC_009495.1     0       56650
+  V       NC_009495.1     8776    8776    1       60      -       tttaaatgaagaaaa NC_009698.1     8776    8791    +
+  V       NC_009495.1     52826   52827   1       60      t       c       NC_009698.1     52841   52842   +
+  V       NC_009495.1     55514   55515   1       60      a       g       NC_009698.1     55529   55530   +
+  V       NC_009495.1     55533   55534   1       60      t       c       NC_009698.1     55548   55549   +
-다음 정보를 정리할 것.
+asm.var.txt에서 R로 시작하는 줄은 하나의 query contig가 덮는 영역을, V로 시작하는 줄은 variant를 뜻한다. paftools.js가 variant calling을 할 때 취하는 alignment의 최소 길이가 정해져 있으므로 short read를 활용할 때에는 조정이 필요하다.
-https://www.biostars.org/p/396962/
-https://github.com/iqbal-lab-org/varifier
-https://genomebiology.biomedcentral.com/articles/10.1186/s13059-021-02473-1
-https://lh3.github.io/2018/04/02/minimap2-and-the-future-of-bwa
-https://academic.oup.com/bioinformatics/article/32/14/2103/1742895