比对软件测试

最后发布时间 : 2023-08-06 19:27:47 浏览量 :

学习资料

完整代码

reads2是reads1的反向互补序列

reads_1.fq

@seq1 match
TGAAGCCAAAGAACAAGATGCGCTAGTGGACAGATTGCTGACCAGGGGCTTGAGAGCTG
+
JIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJJIJJJJIJJIJJJJIJJI

使用下面命令生成reads2

echo "TGAAGCCAAAGAACAAGATGCGCTAGTGGACAGATTGCTGACCAGGGGCTTGAGAGCTG" | rev | tr 'ATCG' 'TAGC'

reads_2.fq

@seq1 match
CAGCTCTCAAGCCCCTGGTCAGCAATCTGTCCACTAGCGCATCTTGTTCTTTGGCTTCA
+
JIJJIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJIJJJJJIJJIJJJJIJJIJJJJI

bowtie2

bowtie2 -x  bowtie2/ref1.fa -1 reads_1.fq -2 reads_2.fq

1 reads; of these:
  1 (100.00%) were paired; of these:
    0 (0.00%) aligned concordantly 0 times
    1 (100.00%) aligned concordantly exactly 1 time
    0 (0.00%) aligned concordantly >1 times
    ----
    0 pairs aligned concordantly 0 times; of these:
      0 (0.00%) aligned discordantly 1 time
    ----
    0 pairs aligned 0 times concordantly or discordantly; of these:
      0 mates make up the pairs; of these:
        0 (0.00%) aligned 0 times
        0 (0.00%) aligned exactly 1 time
        0 (0.00%) aligned >1 times
100.00% overall alignment rate
@HD     VN:1.5  SO:unsorted     GO:query
@SQ     SN:ref1 LN:105
@PG     ID:bowtie2      PN:bowtie2      VN:2.5.1        CL:"/home/wy/workspace/bowtie2-master/bowtie2-align-s --wrapper basic-0 -x bowtie2/ref1.fa -1 reads_1.fq -2 reads_2.fq"
seq1    99      ref1    15      42      59M     =       15      59      TGAAGCCAAAGAACAAGATGCGCTAGTGGACAGATTGCTGACCAGGGGCTTGAGAGCTG     JIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJJIJJJJIJJIJJJJIJJI  AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:59 YS:i:0  YT:Z:CP
seq1    147     ref1    15      42      59M     =       15      -59     TGAAGCCAAAGAACAAGATGCGCTAGTGGACAGATTGCTGACCAGGGGCTTGAGAGCTG     IJJJJIJJIJJJJIJJIJJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJIJ  AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:59 YS:i:0  YT:Z:CP

可以看到bowtie2的flag为

99
- read paired (0x1)
- read mapped in proper pair (0x2)
- mate reverse strand (0x20)
- first in pair (0x40)
147
- read paired (0x1)
- read mapped in proper pair (0x2)
- read reverse strand (0x10)
- second in pair (0x80)

bwa mem

bwa mem bwa/ref1.fa reads_1.fq reads_2.fq

[M::bwa_idx_load_from_disk] read 0 ALT contigs
@SQ     SN:ref1 LN:105
@PG     ID:bwa  PN:bwa  VN:0.7.17-r1188 CL:bwa mem bwa/ref1.fa reads_1.fq reads_2.fq
[M::process] read 2 sequences (118 bp)...
[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (0, 1, 0, 0)
[M::mem_pestat] skip orientation FF as there are not enough pairs
[M::mem_pestat] skip orientation FR as there are not enough pairs
[M::mem_pestat] skip orientation RF as there are not enough pairs
[M::mem_pestat] skip orientation RR as there are not enough pairs
[M::mem_process_seqs] Processed 2 reads in 0.001 CPU sec, 0.001 real sec
seq1    97      ref1    15      60      59M     =       15      59      TGAAGCCAAAGAACAAGATGCGCTAGTGGACAGATTGCTGACCAGGGGCTTGAGAGCTG     JIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJJIJJJJIJJIJJJJIJJI  NM:i:0  MD:Z:59 MC:Z:59M        AS:i:59 XS:i:0
seq1    145     ref1    15      60      59M     =       15      -59     TGAAGCCAAAGAACAAGATGCGCTAGTGGACAGATTGCTGACCAGGGGCTTGAGAGCTG     IJJJJIJJIJJJJIJJIJJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJIJ  NM:i:0  MD:Z:59 MC:Z:59M        AS:i:59 XS:i:0
[main] Version: 0.7.17-r1188
[main] CMD: bwa mem bwa/ref1.fa reads_1.fq reads_2.fq
[main] Real time: 0.002 sec; CPU: 0.005 sec

可以看到bwa mem 的flag为

97
- read paired (0x1)
- mate reverse strand (0x20)
- first in pair (0x40)
145
- read paired (0x1)
- read reverse strand (0x10)
- second in pair (0x80)

reads2与reads1完全相同

reads_1.fq

@seq1 match
TGAAGCCAAAGAACAAGATGCGCTAGTGGACAGATTGCTGACCAGGGGCTTGAGAGCTG
+
JIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJJIJJJJIJJIJJJJIJJI

reads_2_discordantly.fq

@seq1 match
TGAAGCCAAAGAACAAGATGCGCTAGTGGACAGATTGCTGACCAGGGGCTTGAGAGCTG
+
JIJJIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJIJJJJJIJJIJJJJIJJIJJJJI

bowtie2

bowtie2 -x  bowtie2/ref1.fa -1 reads_1.fq -2 reads_2_discordantly.fq

1 reads; of these:
  1 (100.00%) were paired; of these:
    1 (100.00%) aligned concordantly 0 times
    0 (0.00%) aligned concordantly exactly 1 time
    0 (0.00%) aligned concordantly >1 times
    ----
    1 pairs aligned concordantly 0 times; of these:
      0 (0.00%) aligned discordantly 1 time
    ----
    1 pairs aligned 0 times concordantly or discordantly; of these:
      2 mates make up the pairs; of these:
        1 (50.00%) aligned 0 times
        1 (50.00%) aligned exactly 1 time
        0 (0.00%) aligned >1 times
50.00% overall alignment rate
@HD     VN:1.5  SO:unsorted     GO:query
@SQ     SN:ref1 LN:105
@PG     ID:bowtie2      PN:bowtie2      VN:2.5.1        CL:"/home/wy/workspace/bowtie2-master/bowtie2-align-s --wrapper basic-0 -x bowtie2/ref1.fa -1 reads_1.fq -2 reads_2_discordantly.fq"
seq1    73      ref1    15      42      59M     =       15      0       TGAAGCCAAAGAACAAGATGCGCTAGTGGACAGATTGCTGACCAGGGGCTTGAGAGCTG     JIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJJIJJJJIJJIJJJJIJJI  AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:59 YT:Z:UP
seq1    133     ref1    15      0       *       =       15      0       TGAAGCCAAAGAACAAGATGCGCTAGTGGACAGATTGCTGACCAGGGGCTTGAGAGCTG     JIJJIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJIJJJJJIJJIJJJJIJJIJJJJI  YT:Z:UP

可以看到bowtie2的flag为

73
- read paired (0x1)
- mate unmapped (0x8)
- first in pair (0x40)
133
- read paired (0x1)
- read unmapped (0x4)
- second in pair (0x80)

bwa mem

bwa mem bwa/ref1.fa reads_1.fq reads_2_discordantly.fq

[M::bwa_idx_load_from_disk] read 0 ALT contigs
@SQ     SN:ref1 LN:105
@PG     ID:bwa  PN:bwa  VN:0.7.17-r1188 CL:bwa mem bwa/ref1.fa reads_1.fq reads_2_discordantly.fq
[M::process] read 2 sequences (118 bp)...
[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (0, 0, 0, 0)
[M::mem_pestat] skip orientation FF as there are not enough pairs
[M::mem_pestat] skip orientation FR as there are not enough pairs
[M::mem_pestat] skip orientation RF as there are not enough pairs
[M::mem_pestat] skip orientation RR as there are not enough pairs
[M::mem_process_seqs] Processed 2 reads in 0.001 CPU sec, 0.001 real sec
seq1    65      ref1    15      60      59M     =       15      0       TGAAGCCAAAGAACAAGATGCGCTAGTGGACAGATTGCTGACCAGGGGCTTGAGAGCTG     JIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJJIJJJJIJJIJJJJIJJI  NM:i:0  MD:Z:59 MC:Z:59M        AS:i:59 XS:i:0
seq1    129     ref1    15      60      59M     =       15      0       TGAAGCCAAAGAACAAGATGCGCTAGTGGACAGATTGCTGACCAGGGGCTTGAGAGCTG     JIJJIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJIJJJJJIJJIJJJJIJJIJJJJI  NM:i:0  MD:Z:59 MC:Z:59M        AS:i:59 XS:i:0
[main] Version: 0.7.17-r1188
[main] CMD: bwa mem bwa/ref1.fa reads_1.fq reads_2_discordantly.fq
[main] Real time: 0.002 sec; CPU: 0.005 sec

可以看到bwa mem 的flag为

65
- read paired (0x1)
- first in pair (0x40)
129
- read paired (0x1)
- second in pair (0x80)

reads2是reads1的互补序列

reads_1.fq

@seq1 match
TGAAGCCAAAGAACAAGATGCGCTAGTGGACAGATTGCTGACCAGGGGCTTGAGAGCTG
+
JIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJJIJJJJIJJIJJJJIJJI

reads_2_discordantly.fq

@seq1 match
ACTTCGGTTTCTTGTTCTACGCGATCACCTGTCTAACGACTGGTCCCCGAACTCTCGAC
+
JIJJIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJIJJJJJIJJIJJJJIJJIJJJJI

bowtie2

bowtie2 -x  bowtie2/ref1.fa -1 reads_1.fq -2 reads_2_tr.fq

1 reads; of these:
  1 (100.00%) were paired; of these:
    1 (100.00%) aligned concordantly 0 times
    0 (0.00%) aligned concordantly exactly 1 time
    0 (0.00%) aligned concordantly >1 times
    ----
    1 pairs aligned concordantly 0 times; of these:
      0 (0.00%) aligned discordantly 1 time
    ----
    1 pairs aligned 0 times concordantly or discordantly; of these:
      2 mates make up the pairs; of these:
        1 (50.00%) aligned 0 times
        1 (50.00%) aligned exactly 1 time
        0 (0.00%) aligned >1 times
50.00% overall alignment rate
@HD     VN:1.5  SO:unsorted     GO:query
@SQ     SN:ref1 LN:105
@PG     ID:bowtie2      PN:bowtie2      VN:2.5.1        CL:"/home/wy/workspace/bowtie2-master/bowtie2-align-s --wrapper basic-0 -x bowtie2/ref1.fa -1 reads_1.fq -2 reads_2_tr.fq"
seq1    73      ref1    15      42      59M     =       15      0       TGAAGCCAAAGAACAAGATGCGCTAGTGGACAGATTGCTGACCAGGGGCTTGAGAGCTG     JIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJJIJJJJIJJIJJJJIJJI  AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:59 YT:Z:UP
seq1    133     ref1    15      0       *       =       15      0       ACTTCGGTTTCTTGTTCTACGCGATCACCTGTCTAACGACTGGTCCCCGAACTCTCGAC     JIJJIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJIJJJJJIJJIJJJJIJJIJJJJI  YT:Z:UP

可以看到bowtie2的flag为

73
- read paired (0x1)
- mate unmapped (0x8)
- first in pair (0x40)
133
- read paired (0x1)
- read unmapped (0x4)
- second in pair (0x80)

bwa mem

bwa mem bwa/ref1.fa reads_1.fq reads_2_tr.fq

[M::bwa_idx_load_from_disk] read 0 ALT contigs
@SQ     SN:ref1 LN:105
@PG     ID:bwa  PN:bwa  VN:0.7.17-r1188 CL:bwa mem bwa/ref1.fa reads_1.fq reads_2_tr.fq
[M::process] read 2 sequences (118 bp)...
[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (0, 0, 0, 0)
[M::mem_pestat] skip orientation FF as there are not enough pairs
[M::mem_pestat] skip orientation FR as there are not enough pairs
[M::mem_pestat] skip orientation RF as there are not enough pairs
[M::mem_pestat] skip orientation RR as there are not enough pairs
[M::mem_process_seqs] Processed 2 reads in 0.001 CPU sec, 0.001 real sec
seq1    73      ref1    15      60      59M     =       15      0       TGAAGCCAAAGAACAAGATGCGCTAGTGGACAGATTGCTGACCAGGGGCTTGAGAGCTG     JIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJJIJJJJIJJIJJJJIJJI     NM:i:0     MD:Z:59 AS:i:59 XS:i:0
seq1    133     ref1    15      0       *       =       15      0       ACTTCGGTTTCTTGTTCTACGCGATCACCTGTCTAACGACTGGTCCCCGAACTCTCGAC     JIJJIJJJJIJJIJJJJIJJIJJJJIJJIJJJJIJJIJJJJJIJJIJJJJIJJIJJJJI     MC:Z:59M   AS:i:0  XS:i:0
[main] Version: 0.7.17-r1188
[main] CMD: bwa mem bwa/ref1.fa reads_1.fq reads_2_tr.fq
[main] Real time: 0.002 sec; CPU: 0.005 sec

可以看到bwa mem 的flag为

73
- read paired (0x1)
- mate unmapped (0x8)
- first in pair (0x40)
133
- read paired (0x1)
- read unmapped (0x4)
- second in pair (0x80)

可以看到此时bwa mem的flag值与bowtie2相似

BWT算法与后缀树数据格式转换