使用samtools来对sam/bam/cram相互转换

原创

KeepLearningAI 2023-01-04 11:00:01 博主文章分类：Samtools ©著作权

文章标签 bam sam cram hadoop jar 文章分类 软件研发

©著作权归作者所有：来自51CTO博客作者KeepLearningAI的原创作品，请联系作者获取转载授权，否则将追究法律责任

使用samtools来对sam/bam/cram相互转换

1.sam <=>bam

samtools view -h NA12878.bam >NA12878_2.sam  
samtools view -b -S NA12878.sam > NA12878_2.bam

2. cram=>bam

samtools view -bS artificial.sam >artificial.bam
samtools view -bS artificial.sam >artificial.bam
samtools view -bS artificial.cram >artificial2.bam

遇到问题：

hadoop@Mcnode6:~/cloud/adam/xubo/1000genomes/GIH/NA21144/alignment$ samtools view -b -S NA21144.alt_bwamem_GRCh38DH.20150718.GIH.low_coverage.cram >NA21144.alt_bwamem_GRCh38DH.20150718.GIH.low_coverage.bam
[E::cram_populate_ref] mismatching md5sum for downloaded reference.
Failed to populate reference for id 0
Unable to fetch reference #0 9998..119239
Failure to decode slice
[main_samview] truncated file.
hadoop@Mcnode6:~/cloud/adam/xubo/1000genomes/GIH/NA21144/alignment$

hadoop@Mcnode6:~/cloud/adam/xubo/1000genomes/GIH/NA21144/alignment$ ll
total 6991788
drwxr-xr-x 2 hadoop hadoop       4096  3月  9 22:38 ./
drwxr-xr-x 4 hadoop hadoop       4096  3月  8 21:39 ../
-rw-rw-r-- 1 hadoop hadoop     116162  3月  9 22:30 NA21144.alt_bwamem_GRCh38DH.20150718.GIH.low_coverage.bam
-rw-r--r-- 1 hadoop hadoop        877  3月  8 14:19 NA21144.alt_bwamem_GRCh38DH.20150718.GIH.low_coverage.bam.bas
-rw-r--r-- 1 hadoop hadoop 7158968591  3月  8 14:31 NA21144.alt_bwamem_GRCh38DH.20150718.GIH.low_coverage.cram
-rw-r--r-- 1 hadoop hadoop     482019  3月  9 22:38 NA21144.alt_bwamem_GRCh38DH.20150718.GIH.low_coverage.cram.crai

hadoop@Mcnode6:~/cloud/adam/xubo/1000genomes/GIH/NA21144/alignment$ samtools flagstat NA21144.alt_bwamem_GRCh38DH.20150718.GIH.low_coverage2.bam
[W::bam_hdr_read] EOF marker is absent. The input is probably truncated.
248958 + 0 in total (QC-passed reads + QC-failed reads)
481 + 0 secondary
0 + 0 supplementary
4660 + 0 duplicates
247408 + 0 mapped (99.38% : N/A)
248477 + 0 paired in sequencing
124288 + 0 read1
124189 + 0 read2
227967 + 0 properly paired (91.75% : N/A)
245377 + 0 with itself and mate mapped
1550 + 0 singletons (0.62% : N/A)
3892 + 0 with mate mapped to a different chr
1109 + 0 with mate mapped to a different chr (mapQ>=5)
hadoop@Mcnode6:~/cloud/adam/xubo/1000genomes/GIH/NA21144/alignment$ samtools flagstat NA21144.alt_bwamem_GRCh38DH.20150718.GIH.low_coverage3.bam
97098407 + 0 in total (QC-passed reads + QC-failed reads)
179635 + 0 secondary
0 + 0 supplementary
2634031 + 0 duplicates
96638779 + 0 mapped (99.53% : N/A)
96918772 + 0 paired in sequencing
48457840 + 0 read1
48460932 + 0 read2
93116714 + 0 properly paired (96.08% : N/A)
95999516 + 0 with itself and mate mapped
459628 + 0 singletons (0.47% : N/A)
1495404 + 0 with mate mapped to a different chr
565190 + 0 with mate mapped to a different chr (mapQ>=5)

3.bam<=>cram:

samtools view -C -T ref.fa aln.bam > aln.cram

java -jar cramtools-3.0.jar bam -O yeast.bam -I yeast.cram -R yeast.fasta

java -jar cramtools-3.0.jar cram -O yeast2.cram  -I yeast.bam -R yeast.fasta

运行记录：

hadoop@Mcnode6:~/cloud/adam/xubo/yeast201603101125/test1$ java -jar cramtools-3.0.jar bam -O yeast.bam -I yeast.cram -R yeast.fasta 
hadoop@Mcnode6:~/cloud/adam/xubo/yeast201603101125/test1$ ls
cramtools-3.0.jar  yeast.bam  yeast.cram  yeast.fasta
hadoop@Mcnode6:~/cloud/adam/xubo/yeast201603101125/test1$ ll
total 19008
drwxrwxr-x 2 hadoop hadoop     4096  3月 10 15:23 ./
drwxrwxr-x 3 hadoop hadoop     4096  3月 10 15:01 ../
-rw-rw-r-- 1 hadoop hadoop  3986091  3月 10 15:01 cramtools-3.0.jar
-rw-rw-r-- 1 hadoop hadoop  2130246  3月 10 15:25 yeast.bam
-rw-rw-r-- 1 hadoop hadoop   967382  3月 10 15:01 yeast.cram
-rw-rw-r-- 1 hadoop hadoop 12360755  3月 10 15:23 yeast.fasta
hadoop@Mcnode6:~/cloud/adam/xubo/yeast201603101125/test1$ java -jar cramtools-3.0.jar bam -O yeast2.bam -I yeast.cram -R yeast.fasta 
hadoop@Mcnode6:~/cloud/adam/xubo/yeast201603101125/test1$ ls
cramtools-3.0.jar  yeast2.bam  yeast.bam  yeast.cram  yeast.fasta
hadoop@Mcnode6:~/cloud/adam/xubo/yeast201603101125/test1$ ll
total 21092
drwxrwxr-x 2 hadoop hadoop     4096  3月 10 15:26 ./
drwxrwxr-x 3 hadoop hadoop     4096  3月 10 15:01 ../
-rw-rw-r-- 1 hadoop hadoop  3986091  3月 10 15:01 cramtools-3.0.jar
-rw-rw-r-- 1 hadoop hadoop  2130242  3月 10 15:26 yeast2.bam
-rw-rw-r-- 1 hadoop hadoop  2130246  3月 10 15:25 yeast.bam
-rw-rw-r-- 1 hadoop hadoop   967382  3月 10 15:01 yeast.cram
-rw-rw-r-- 1 hadoop hadoop 12360755  3月 10 15:23 yeast.fasta
hadoop@Mcnode6:~/cloud/adam/xubo/yeast201603101125/test1$ java -jar cramtools-3.0.jar bam -O yeast3.bam -I yeast.cram -R yeast.fasta 
hadoop@Mcnode6:~/cloud/adam/xubo/yeast201603101125/test1$ ll
total 23176
drwxrwxr-x 2 hadoop hadoop     4096  3月 10 15:28 ./
drwxrwxr-x 3 hadoop hadoop     4096  3月 10 15:01 ../
-rw-rw-r-- 1 hadoop hadoop  3986091  3月 10 15:01 cramtools-3.0.jar
-rw-rw-r-- 1 hadoop hadoop  2130242  3月 10 15:26 yeast2.bam
-rw-rw-r-- 1 hadoop hadoop  2130242  3月 10 15:28 yeast3.bam
-rw-rw-r-- 1 hadoop hadoop  2130246  3月 10 15:25 yeast.bam
-rw-rw-r-- 1 hadoop hadoop   967382  3月 10 15:01 yeast.cram
-rw-rw-r-- 1 hadoop hadoop 12360755  3月 10 15:23 yeast.fasta
hadoop@Mcnode6:~/cloud/adam/xubo/yeast201603101125/test1$ java -jar cramtools-3.0.jar bam -O yeast4.bam -I yeast.cram -R yeast.fasta 
hadoop@Mcnode6:~/cloud/adam/xubo/yeast201603101125/test1$ ll
total 25260
drwxrwxr-x 2 hadoop hadoop     4096  3月 10 15:30 ./
drwxrwxr-x 3 hadoop hadoop     4096  3月 10 15:01 ../
-rw-rw-r-- 1 hadoop hadoop  3986091  3月 10 15:01 cramtools-3.0.jar
-rw-rw-r-- 1 hadoop hadoop  2130242  3月 10 15:26 yeast2.bam
-rw-rw-r-- 1 hadoop hadoop  2130242  3月 10 15:28 yeast3.bam
-rw-rw-r-- 1 hadoop hadoop  2130242  3月 10 15:31 yeast4.bam
-rw-rw-r-- 1 hadoop hadoop  2130246  3月 10 15:25 yeast.bam
-rw-rw-r-- 1 hadoop hadoop   967382  3月 10 15:01 yeast.cram
-rw-rw-r-- 1 hadoop hadoop 12360755  3月 10 15:23 yeast.fasta
hadoop@Mcnode6:~/cloud/adam/xubo/yeast201603101125/test1$ java -jar cramtools-3.0.jar 
Version 3.0-b48

Usage: cramtools [options] [command] [command options]
  Options:    -h, --help  Print help and quit (default: false)
  Commands:
    bam         CRAM to BAM conversion. 
    cram        BAM to CRAM converter. 
    index       BAM/CRAM indexer. 
    merge       Tool to merge CRAM or BAM files. 
    fastq       CRAM to FastQ dump conversion. 
    fixheader   A tool to fix CRAM header without re-writing the whole file.
    getref      Download reference sequences.
    qstat       Quality score statistics.

hadoop@Mcnode6:~/cloud/adam/xubo/yeast201603101125/test1$ java -jar cramtools-3.0.jar cram
Version 3.0-b48

Usage: <main class> [options]
 
  Options:    --capture-all-tags              Capture all tags. (default: false)
    --capture-tags                  Capture the tags listed, for example 'OQ:XA:XB' (default: )
    --encrypt                       Encrypt the CRAM file. (default: false)
    --ignore-md5-mismatch           Fail on MD5 mismatch if true, or correct (overwrite) the checksums and continue if false. (default: false)
    --ignore-tags                   Ignore the tags listed, for example 'OQ:XA:XB' (default: )
    --inject-sq-uri                 Inject or change the @SQ:UR header fields to point to ENA reference service.  (default: false)
    --input-bam-file, -I            Path to a BAM file to be converted to CRAM. Omit if standard input (pipe).
    --input-is-sam                  Input is in SAM format. (default: false)
    --lossless-quality-score, -Q    Preserve all quality scores. Overwrites '--lossless-quality-score'. (default: false)
    --lossy-quality-score-spec, -L  A string specifying what quality scores should be preserved. (default: )
    --max-records                   Stop after compressing this many records.  (default: 9223372036854775807)
    --output-cram-file, -O          The path for the output CRAM file. Omit if standard output (pipe).
    --preserve-read-names, -n       Preserve all read names. (default: false)
    --reference-fasta-file, -R      The reference fasta file, uncompressed and indexed (.fai file, use 'samtools faidx'). 
    -h, --help                      Print help and quit (default: false)
    -l, --log-level                 Change log level: DEBUG, INFO, WARNING, ERROR. (default: ERROR)

使用上述指令转换后明显有数据压缩：

hadoop@Mcnode6:~/cloud/adam/xubo/yeast201603101125/test1$ ll
total 26760
drwxrwxr-x 2 hadoop hadoop     4096  3月 10 16:49 ./
drwxrwxr-x 3 hadoop hadoop     4096  3月 10 15:01 ../
-rw-rw-r-- 1 hadoop hadoop  3986091  3月 10 15:01 cramtools-3.0.jar
-rw-rw-r-- 1 hadoop hadoop  2130242  3月 10 15:26 yeast2.bam
-rw-rw-r-- 1 hadoop hadoop   510298  3月 10 15:33 yeast2.cram
-rw-rw-r-- 1 hadoop hadoop  2130242  3月 10 15:28 yeast3.bam
-rw-rw-r-- 1 hadoop hadoop   510301  3月 10 15:40 yeast3.cram
-rw-rw-r-- 1 hadoop hadoop  2130242  3月 10 15:31 yeast4.bam
-rw-rw-r-- 1 hadoop hadoop   510301  3月 10 16:50 yeast5.cram
-rw-rw-r-- 1 hadoop hadoop  2130246  3月 10 15:25 yeast.bam
-rw-rw-r-- 1 hadoop hadoop   967382  3月 10 15:01 yeast.cram
-rw-rw-r-- 1 hadoop hadoop 12360755  3月 10 15:23 yeast.fasta

hadoop@Mcnode6:~/cloud/adam/xubo/yeast201603101125/test1$ samtools view yeast5.cram |head -20
SRR507778.19213 147 I 62  60  36M = 3183  3087  ATCCTAACACTACCCTAACACAGCCCTAATCTAACC  * MD:Z:36 NM:i:0
SRR507778.12312 147 I 205 60  36M = 3626  3387  CCACTCACCCACCGTTACCCTCCAATTACCCATATC  * MD:Z:36 NM:i:0
SRR507778.11604 83  I 402 60  36M = 3869  3503  CTCACTTGTATACTGATTTTACGTACGCACACGGAT  * MD:Z:36 NM:i:0
SRR507778.10609 83  I 2661  40  36M = 6131  3506  TGAATTCGTACAACATTAAACGTGTGTTGGGAGTCG  * MD:Z:36 NM:i:0
SRR507778.6249  147 I 2925  60  36M = 6404  3445  TTTCTAAGTGGGATTTTTCTTAATCCTTGGATTCTT  * MD:Z:36 NM:i:0
SRR507778.14609 129 I 3048  60  36M IV  1525643 0 AAAAGTAGCCGTTCATTTCCCTTCCGATTTCATTCC  * MD:Z:36 NM:i:0
SRR507778.20233 83  I 3132  60  36M = 6388  3292  TATTTGTGTCCCATTCTCGTAGATAAAATTCTTGGA  * MD:Z:36 NM:i:0
SRR507778.19213 99  I 3183  60  36M = 62  -3087 ATTTTCTTCATAAAGAAGCTTTCAAGATATAAGATA  * MD:Z:36 NM:i:0
SRR507778.20882 73  I 3259  60  36M = 3259  0 CAAAAAGGAAAGCATGGAGGGAAACAGTAAACAGTG  * MD:Z:36 NM:i:0
SRR507778.20882 133 I 3259  0 * = 3259  0 GTGGTGTGTGTGGGTGAGGTGTGGGTGTGGGGAGGG  *
SRR507778.12312 99  I 3626  60  36M = 205 -3387 GTATCTGATGTTTTTTTAGTAATTTCTTTGTAAATA  * MD:Z:36 NM:i:0
SRR507778.11604 163 I 3869  60  36M = 402 -3503 TTTTTGAAAATATTCTGAGGTAAAAGCCATTAAGGT  * MD:Z:36 NM:i:0
SRR507778.24515 83  I 4004  60  36M = 7814  3846  GATGTTTCAAGGCCTGAAGTTTGAATATTTATGTAG  * MD:Z:36 NM:i:0
SRR507778.19471 83  I 4627  60  36M = 8153  3562  GGCAGAGTTTCCAAAAAAAATTGTTAATCGACAAAG  * MD:Z:36 NM:i:0
SRR507778.15626 83  I 4748  60  36M = 8861  4149  TTTAAATTGTATTGAGTGCTTCAGTCATTGCAAAAT  * MD:Z:36 NM:i:0
SRR507778.7265  147 I 4894  60  36M = 8228  3300  TATCTATCACAAAGGAGACAAAATCGTTGATAAAAA  * MD:Z:36 NM:i:0
SRR507778.14364 83  I 5516  60  36M = 9133  3653  TATGATATAAAAACTCGGACCCTGTTTTACTTCTTT  * MD:Z:36 NM:i:0
SRR507778.10609 163 I 6131  60  36M = 2661  -3506 CATACGTTGATTAGTACTGTTGGTCTCTCATTGAAA  * MD:Z:36 NM:i:0
SRR507778.20233 163 I 6388  60  36M = 3132  -3292 ACCAATTTGACGTTAATTTTAAATGCGTTCTGAAGT  * MD:Z:36 NM:i:0
SRR507778.6249  99  I 6404  60  36M = 2925  -3445 TTTTAAATGCGTTCTGAAGTTTCTTAAATAACCCGG  * MD:Z:36 NM:i:0
hadoop@Mcnode6:~/cloud/adam/xubo/yeast201603101125/test1$ samtools view yeast.cram |head -20
SRR507778.19213 147 I 62  60  36M = 3183  3087  ATCCTAACACTACCCTAACACAGCCCTAATCTAACC  15=@9:@C3<CBGGGDGDBGDFCC?>GGG<GGGDGG AS:i:36 XS:i:19 MD:Z:36 NM:i:0
SRR507778.12312 147 I 205 60  36M = 3626  3387  CCACTCACCCACCGTTACCCTCCAATTACCCATATC  GD<B?B>DGGBGGGGGCIIIIGIIIIIEIIGIIIII  AS:i:36 XS:i:26 MD:Z:36 NM:i:0
SRR507778.11604 83  I 402 60  36M = 3869  3433  CTCACTTGTATACTGATTTTACGTACGCACACGGAT  GHGHIIIIIIIIIIIHIIIIIIIIIIIIIIIIIIII  AS:i:36 XS:i:0  MD:Z:36 NM:i:0
SRR507778.10609 83  I 2661  40  36M = 6131  3436  TGAATTCGTACAACATTAAACGTGTGTTGGGAGTCG  IIIGFIIGIIHIIIIGIIIIIIIIIIHIIIIIIGII  AS:i:36 XS:i:36 MD:Z:36 NM:i:0
SRR507778.6249  147 I 2925  60  36M = 6404  3445  TTTCTAAGTGGGATTTTTCTTAATCCTTGGATTCTT  GGGGADIGIHIIHHHIEHIHIHI<DIIIIIIIGIIF AS:i:36 XS:i:0  MD:Z:36 NM:i:0
SRR507778.14609 129 I 3048  60  36M IV  1525643 0 AAAAGTAGCCGTTCATTTCCCTTCCGATTTCATTCC  >5833+?=8>B@FBF?9B7AGGGB<G@BGGGGEGD>  AS:i:36 XS:i:0  MD:Z:36 NM:i:0
SRR507778.20233 83  I 3132  60  36M = 6388  3222  TATTTGTGTCCCATTCTCGTAGATAAAATTCTTGGA  IIIIIIIGGIIIIIIIIIIIIIIIIIIIIIIBIIII  AS:i:36 XS:i:0  MD:Z:36 NM:i:0
SRR507778.19213 99  I 3183  60  36M = 62  -3087 ATTTTCTTCATAAAGAAGCTTTCAAGATATAAGATA  HHHGHGDAHHHHHEHHHHHHHHGHEGBDGGGGG<GE AS:i:36 XS:i:0  MD:Z:36 NM:i:0
SRR507778.20882 73  I 3259  60  36M = 3259  0 CAAAAAGGAAAGCATGGAGGGAAACAGTAAACAGTG  @GGGGGG>GGBD4DDGGEDGDDG@GAA1CBEEEE3D AS:i:36 XS:i:0  MD:Z:36 NM:i:0
SRR507778.20882 133 I 3259  0 * = 3259  0 GTGGTGTGTGTGGGTGAGGTGTGGGTGTGGGGAGGG  EGBG8GCB8BBBB#######################  AS:i:0  XS:i:0
SRR507778.12312 99  I 3626  60  36M = 205 -3387 GTATCTGATGTTTTTTTAGTAATTTCTTTGTAAATA  IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIHHI  AS:i:36 XS:i:0  MD:Z:36 NM:i:0
SRR507778.11604 163 I 3869  60  36M = 402 -3433 TTTTTGAAAATATTCTGAGGTAAAAGCCATTAAGGT  IIIIIIHIIIIIIIIIIIIIIIIIIIIHIIIIHIIE  AS:i:36 XS:i:0  MD:Z:36 NM:i:0
SRR507778.24515 83  I 4004  60  36M = 7814  3776  GATGTTTCAAGGCCTGAAGTTTGAATATTTATGTAG  IHIIIIIIIIIIIIHIIIIHIIIIIIIIIIIIIIII  AS:i:36 XS:i:0  MD:Z:36 NM:i:0
SRR507778.19471 83  I 4627  60  36M = 8153  3492  GGCAGAGTTTCCAAAAAAAATTGTTAATCGACAAAG  HGHBHGGHGHHHGHHHHHHHHGGGGGDDD=BDGGGG  AS:i:36 XS:i:20 MD:Z:36 NM:i:0
SRR507778.15626 83  I 4748  60  36M = 8861  4079  TTTAAATTGTATTGAGTGCTTCAGTCATTGCAAAAT  IIIIIIIIIIIIIIIIIIIIIIHIIIIIIIIIIIII  AS:i:36 XS:i:0  MD:Z:36 NM:i:0
SRR507778.7265  147 I 4894  60  36M = 8228  3300  TATCTATCACAAAGGAGACAAAATCGTTGATAAAAA  GGBDGGIIIFIGHHIIDDGGGGGDGDIIDIIEIIII  AS:i:36 XS:i:0  MD:Z:36 NM:i:0
SRR507778.14364 83  I 5516  60  36M = 9133  3583  TATGATATAAAAACTCGGACCCTGTTTTACTTCTTT  IIIIIIIHIIIIIIIIIIIIIIIIIIIIIIIIIIII  AS:i:36 XS:i:0  MD:Z:36 NM:i:0
SRR507778.10609 163 I 6131  60  36M = 2661  -3436 CATACGTTGATTAGTACTGTTGGTCTCTCATTGAAA  HIHIHIGIIIHIIIHIHIIIGIIIGIEHDIIIHIHG  AS:i:36 XS:i:0  MD:Z:36 NM:i:0
SRR507778.20233 163 I 6388  60  36M = 3132  -3222 ACCAATTTGACGTTAATTTTAAATGCGTTCTGAAGT  IIIIIIIIIIIIIIIIIIIIIIIIIIIIIGIIIHII  AS:i:36 XS:i:0  MD:Z:36 NM:i:0
SRR507778.6249  99  I 6404  60  36M = 2925  -3445 TTTTAAATGCGTTCTGAAGTTTCTTAAATAACCCGG  GGGGGGGGGGHHHHHHGDHHHHGHHFGHHHGFGGGG  AS:i:36 XS:i:0  MD:Z:36 NM:i:0
hadoop@Mcnode6:~/cloud/adam/xubo/yeast201603101125/test1$