将微基因全基因组cram文件解压转换为bam文件教程

转成bam的必要性是啥呢

不好意思，我在转换的时候可能对参考基因组数据理解有误解，感谢楼上同学为我指出错误。大家可以直接用楼上同学的操作方法去解压cram文件。

楼主你好像理解错了-T参数的意思。
 
-T 参数其实是需要一个参考基因组ref文件，这个在cram/bam的header里能看到（samtools view -H input.cram | less ），WeGene用的应该是这个版本（human_g1k_v37.fasta），去下载一个就好了。
 
1. 下载ref文件
wget http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/human_g1k_v37.fasta.gz
2. 解压
gunzip human_g1k_v37.fasta.gz
3. 给ref建个index
samtools faidx human_g1k_v37.fasta
4. 把cram转成bam就好了
samtools view -b -T human_g1k_v37.fasta -o output.bam input.cram

建议通过管道命令行运行，不然中间文件太大了

支持

有大佬来了，赞一个以支持一下

如果有用请点赞支持，申请加精！

8.接下来把bam和bai文件放到一个目录下，用IGV_2.16.0可以打开你的基因序列
 
 


 
补充一些关于SAM和BAM文件的基本知识，这可能有助于更好地理解如何将BAM文件转换为SAM文件以及它们之间的区别：
SAM（Sequence Alignment/Map）是一种文本格式，用于描述序列比对的结果，包括比对到参考序列的序列读取、比对位置、质量值等信息。SAM 文件中的每一行表示一个比对结果，以制表符分隔的字段包括序列名称、比对标志、参考序列名称、比对起始位置、质量值、CIGAR字符串等。SAM 格式是一种人类可读的格式，方便进行手动编辑或解析，但是文件体积比较大，不适合存储大规模的测序数据。
 

7.最后我们需要生成bai引索文件来加速打开bam文件的速度
分别执行命令：
 
samtools index -@ 8 63288812349876.Y.bam
 
samtools index -@ 8 63288812349876.bam
 
上面的命令是生成Y染色体的bai文件
下面的命令是生成cram文件的bai文件
到现在你就获得了原始的bam测序文件，我的文件大小63G，你们的多大我就不清楚了。
 
 

6.获得fastA文件后，我们就可以用fastA文件去作为校对列，去解压cram文件，将cram文件转换为bam文件，我本人做的是30X的全基因组序列，预测结果为900G文件大小。请提前准备好对应大小的移动硬盘来存储数据。
 
需要注意的是，微基因的参考基因组序列是GRCh37/hg19,请自行下载对应的参考序列，使用不同的参考序列会导致结果不准确。如果不知道如何下载可以请教地表最强AI，ChatGPT。

这里修改一下，经过热心网友的讨论，我们可以直接使用现成的fasta文件而不必自己提取fasta
1. 下载ref文件
 
wget http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/human_g1k_v37.fasta.gz
 
2. 解压
我直接在windows11下使用360压缩进行了解压，因为是虚拟机共享文件夹，所以省去了命令指令

3. 给ref建个index
samtools faidx human_g1k_v37.fasta
 
4. 把cram转成bam就好了

执行命令:
samtools view -@ 20 -b -T human_g1k_v37.fasta -o 63288812349876.bam 63288812349876.cram
 
 
本命令使用20线程加速转换，如果你的设备很强，请自行修改线程数量。我的32G内存已经全部占满了，CPU占用倒不是很大，只有32%
 






（这个命令会从输入的 CRAM 文件 63288812349876.cram 中解压缩并转换成 BAM 格式的文件 63288812349876.bam。参考基因组文件 63288812349876.Y.bam.sam.1.fasta 用于提供基因组信息，确保输出 BAM 文件中的 reads 在正确的位置上。这里的 -T 选项指定参考基因组文件，-b 选项将输出文件转换为 BAM 格式。
 
请注意，如果使用 CRAM 文件，需要提供参考基因组文件以将 CRAM 文件解压缩成 BAM 文件。但如果使用 BAM 文件，则不需要提供参考基因组文件，因为 BAM 文件已经包含了基因组信息。）
 
 
等到几小时后，你就会获得一个63288812349876.bam的文件，这个文件就是cram解压后的原始数据，到这里你就获得了你要的数据。
我转换花了45分钟，加压后文件大小为62.3GB。

5.获得fastQ文件后,如果您想要使用 fastQ 文件作为参考基因组文件，则需要先将 fastQ 文件转换为 fastA 格式的文件，然后再将其传递给 -T 参数。
samtools不支持将FASTQ文件直接转换为FASTA文件，因为这是两种不同的文件格式，FASTQ文件包含了序列的质量信息，而FASTA文件只包含序列信息。因此，需要使用其他工具来进行转换。
可以使用seqtk工具将FASTQ文件转换为FASTA文件，并且seqtk支持多线程操作。
其中，-a参数表示转换为FASTA格式，-j参数指定使用的线程数，input.fastq为输入的FASTQ文件，output.fasta为输出的FASTA文件。
 
首先安装必要命令库
 
执行命令：
 
安装 zlib 库：
sudo apt-get install zlib1g-dev
 
安装 GCC 编译器：
sudo apt-get install gcc
 
安装seqtk：
git clone https://github.com/lh3/seqtk.git;
cd seqtk; make
 
出现这个提示说明编译 Seqtk 成功了，并生成了一个名为 "seqtk" 的可执行文件。

执行命令：
seqtk seq -a 63288812349876.Y.bam.sam.1.fastq > 63288812349876.Y.bam.sam.1.fasta
 
 
 
该命令将fastQ转换为fastA
执行结束后将得到一个新文件：
63288812349876.Y.bam.sam.1.fasta
 

 

需要等待几分钟，计算机返回执行结果
[M::bam2fq_mainloop] discarded 0 singletons
[M::bam2fq_mainloop] processed 5090624 reads
 


（这是samtools bam2fq命令输出的一部分信息。其中，“[M::bam2fq_mainloop] discarded 0 singletons”表示在转换过程中没有丢弃任何单独的reads。
而“[M::bam2fq_mainloop] processed 5090624 reads”则表示该命令已经处理了5090624个reads。
这个信息告诉我们samtools bam2fq已经成功地将BAM文件转换为了FASTQ文件，并且输出了多少reads。如果你需要的话，你可以根据这个信息检查转换结果是否正确，并在处理大规模数据时，通过这个信息来了解转换进度。）
 
然后获得2个新文件
63288812349876.Y.bam.sam.1.fastq
63288812349876.Y.bam.sam.2.fastq
 

4.然后我们需要从63288812349876.Y.bam.sam文件提取fastQ文件（fastQ是一种存储了生物序列（通常是核酸序列）以及相应的质量评价的文本格式）。
 
通常情况下，从SAM文件中提取fastQ格式的序列数据时，会生成两个fastQ文件，一个包含读1的序列数据，另一个包含读2的序列数据。这是因为Illumina测序数据通常是成对的，即包含了从同一个基因片段的两端读取的序列数据。因此，您需要分别提取这两端的序列，并将其保存在不同的文件中。
 
 
执行命令：
samtools fastq -@ 8 63288812349876.Y.bam.sam -1 63288812349876.Y.bam.sam.1.fastq -2 63288812349876.Y.bam.sam.2.fastq
 
本命令使用8线程加速操作，
 

我现在还差两个WGS没上传NCBI就是这个东西搞的。 很烦

我在家试过的，但是家里的机器好像有点问题，每次快要结束的时候就会弹出报错。我觉得微基因应该提供一下

3.解压cram文件需要使用fastA文件作为参考文件，
提取fastA文件的方法是从63288812349876.Y.bam中提取fastQ,再将fastQ转换为fastA
而bam文件包含fastQ文件，所以
首先需要将63288812349876.Y.bam文件转换为sam文件
 
执行命令
samtools view -@ 8 -h 63288812349876.Y.bam > 63288812349876.Y.bam.sam
将bam文件转换为sam文件，其中使用8线程加速操作
 
然后获得1个新文件
63288812349876.Y.bam.sam