在Windows系统中使用BWA更换全基因组数据参考的大致方法

部分指令参考了WGS Extract V3临时生成的shell脚本。接下来这里直接展示步骤：

一、下载工具与下载参考基因组fasta

如果您使用的Windows系统至少为Windows 10版本，那么进入WGS Extract V3网站，下载Win10tools。

WGS Extract | WGSExtract.github.ioWGS Extract WWW homehttps://wgsextract.github.io/Win10tools工具使用方法以及参考基因组的激活使用方法在下文中也有介绍，用Cygwin窗口或Windows终端执行，本文不再赘述：

【教程】如何使用BCFtools提取全基因组数据到芯片模拟数据？_风澪宙的博客-CSDN博客_bcftools如果您使用的是Windows系统（至少Windows 10），操作这部分的内容建议下载解压Cygwin衍生的win10tools，并打开Cygwin.bat来在Windows模拟部分Linux系统环境。文末有win10tools的下载链接，如果失效，那么需从以下链接下翻找到Win10 Release Users来下载win10tools：WGS Extract Version 3 Beta | WGSExtract.github.iohttps://wgsextract.github.io/ 接着打开使用https://blog.csdn.net/qq_53947118/article/details/125542406

win10tools的网盘下载链接https://pan.baidu.com/s/16f7YU2e9vNOOpHuPvoQXKg?pwd=cyg1

二、观察所用参考基因组

1. 读取全基因组数据文件的头部

samtools view -H --no-PG xxx.bam

或

samtools view -H --no-PG xxx.bam >xxx.txt

2. 直接从终端窗口或打开生成的文本文档，查看自己的数据所用参考基因组的类型（比如图中的为Human_G1K_V37）。

三、转化

前提：

1. 数十GB的全基因组数据在转化时需要消耗大量的硬盘空间，建议留出几百GB甚至1TB以上的空余空间。

2. 转化全基因组数据坐标所花费的时间极长，建议在空闲时间且温度较凉的环境下操作（比如至少连续2天）。

3. 作测试的案例使用了包含的hs38DH的参考基因组，非特殊情况下不建议使用该参考来转。参考基因组的选择要经过谨慎思考，数据质量高的为极佳选择，以下链接仅展示hs37d5、hs38、hs38d1、hs38s（hs38d1+22_KI270879v1_alt，Sequencing.com所用参考）参考基因组的文件下载链接，建议根据需要选择。

hs37d5https://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gzhs38https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gzhs38d1https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids/GCA_000001405.15_GRCh38_no_alt_plus_hs38d1_analysis_set.fna.gzhs38shttps://api.onedrive.com/v1.0/shares/s!AgorjTSMFYpjgR0QualUlHx53-0U/root/content4. 接下来的步骤使用被切裁过的部分全基因组文件来展示，因此文件大小远小于完整的全序文件，实际上这些步骤对完整的全基因组数据也有效。

步骤：

1. 为参考基因组建立bwa专用索引（“xx”可替换为自己所用的参考名）

(1) （可选）为参考基因组建立字典

samtools dict xx.fa.gz > xx.dict

(2) 解压（因为原参考可能被gzip压缩而不兼容其他工具），然后用bgzip重新压缩

bgzip -d xx.fa.gz
bgzip -f xx.fa

(3) 用samtools和bwa建立索引（后者需要等待较长的时间）

samtools index xx.fa.gz
bwa index xx.fa.gz

2. bam/cram文件转化为fastq（全序数据以文件名“Filename00000”为例，且参考以hs38DH为例）

samtools view -uh --no-PG Filename00000.bam |\
samtools sort -n -T "./temp/" -m 1000M -@ 16 -O sam |\
samtools fastq -1 Filename00000_R1.fastq.gz -2 Filename00000_R2.fastq.gz -0 /dev/null -s Filename00000.fastq.gz -n -@ 16

或

samtools view -uh --no-PG -T "./Reference/human_g1k_v37.fasta.gz" Filename00000.cram |\
samtools sort -n -T "./temp/" -m 1000M -@ 16 -O sam |\
samtools fastq -1 Filename00000_R1.fastq.gz -2 Filename00000_R2.fastq.gz -0 /dev/null -s Filename00000.fastq.gz -n -@ 16

3. 使用bwa工具remap全基因组数据（直接生成的话会产生占用数百GB空间的sam文件，因此此处结合bwa mem和samtools fixmate指令）

(1) remap配对的序列

bwa mem -t 16 -B 4 -M -R "@RG\tID:xxxxx\tSM:Filename00000\tPL:xxxxx\tLB:xxxxx" "./Reference/GRCh38_full_analysis_set_plus_decoy_hla.fa" Filename00000_R1.fastq.gz Filename00000_R2.fastq.gz |\
samtools fixmate -m -O bam -@ 16 - Filename00000-hs38DH_PE_fixmate.bam

(2) remap未配对的序列（若未配对数据部分为0B大小，也可省略这一步）

bwa mem -t 16 -B 4 -M -R "@RG\tID:xxxxx\tSM:Filename00000\tPL:xxxxx\tLB:xxxxx" "./Reference/GRCh38_full_analysis_set_plus_decoy_hla.fa" Filename00000.fastq.gz |\
samtools fixmate -m -O bam -@ 16 - Filename00000-hs38DH_SE_fixmate.bam

(3) 合并原始bam文件（若未配对数据部分为0B大小，也可省略这一步）

samtools merge -@ 16 -r Filename00000-hs38DH_fixmate.bam Filename00000-hs38DH_PE_fixmate.bam Filename00000-hs38DH_SE_fixmate.bam

4. 优化bam文件，并建立索引

(1) 排序

samtools sort -T "./temp/" -@ 16 -m 1000M Filename00000-hs38DH_fixmate.bam > Filename00000-hs38DH_sort.bam

(2) 标记重复的reads并输出最终结果（感兴趣的话也可使用dedup去除重复的reads）

samtools markdup -@ 16 -O bam Filename00000-hs38DH_sort.bam - | samtools view -bh -@ 16 - -o Filename00000-hs38DH.bam

或

samtools markdup -@ 16 -O bam Filename00000-hs38DH_sort.bam - | samtools view -Ch -T "./Reference/GRCh38_full_analysis_set_plus_decoy_hla.fa" -@ 16 - > Filename00000-hs38DH.cram

(2) 建立索引

samtools index Filename00000-hs38DH.bam

或

samtools index Filename00000-hs38DH.cram

5. （可选）读取转化好的bam/cram文件信息

(1) 读取文件头部（同上）

samtools view -H --no-PG xxx.bam >xxx_header.txt

或

samtools view -H --no-PG xxx.cram >xxx_header.txt

(2) 计算reads读取情况

samtools idxstats -@ 16 xxx.bam > xxx_idxstats.txt

或

samtools idxstats -@ 16 xxx.cram > xxx_idxstats.txt

(3) 计算数据覆盖程度

samtools coverage xxx.bam > xxx_coverage.txt

或

samtools coverage --reference xx.fa.gz xxx.cram > xxx_coverage.txt

拓展阅读：
FASTQ to BAM / CRAMSamtoolshttps://www.htslib.org/workflow/fastq.htmlWhich human reference genome to use?https://lh3.github.io/2017/11/13/which-human-reference-genome-to-useIt's finally finished! – Genome Informatics SectionToday is a big day. One that was 30+ years in the making. We have finally uncovered every last bit of the human genome! I wanted to celebrate by writing a behind the scenes tribute to the amazing members of the T2T consortium, but that will have to wait. The past two years have been a whirlwind and I need some time to recuperate. For now, just the basics: links to the assembly, the browser, and the papers. And don’t worry, we didn’t forget chrY this time!https://genomeinformatics.github.io/CHM13v2/?msclkid=932a6a2cb29611ecaab4de173e3c0d88