NA12878样本下载和处理
NA12878样本常被用来作为benchmark的参考,介绍一下该样本基因组序列的下载和处理。
1. 下载
- 下载NA12878的cram文件(该网站上只有cram文件,需要后续根据使用的参考基因组转成bam文件):https://www.internationalgenome.org/data-portal/sample/NA12878
#!/bin/bash
# all file download from https://www.internationalgenome.org/data-portal/sample/NA12878
## [1] low coverage WGS data
# cram and crai
wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000_genomes_project/data/CEU/NA12878/alignment/NA12878.alt_bwamem_GRCh38DH.20150718.CEU.low_coverage.cram
# fasta file
# wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR622/SRR622461/SRR622461_2.fastq.gz
# wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR622/SRR622461/SRR622461_1.fastq.gz
- 下载参考基因组(一定要下载对应的,否则转成的bam文件在后续运行时会出现文件不完整的报错)
wget -c https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/20150713_location_of_centromeres_and_other_regions.txt /home/data/sda/tzy/NA12878/GRCh38_reference_genome
wget -c https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla-extra.fa /home/data/sda/tzy/NA12878/GRCh38_reference_genome
wget -c https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.dict /home/data/sda/tzy/NA12878/GRCh38_reference_genome
wget -c https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa /home/data/sda/tzy/NA12878/GRCh38_reference_genome
wget -c https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa.alt /home/data/sda/tzy/NA12878/GRCh38_reference_genome
wget -c https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa.amb /home/data/sda/tzy/NA12878/GRCh38_reference_genome
wget -c https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa.ann /home/data/sda/tzy/NA12878/GRCh38_reference_genome
wget -c https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa.bwt /home/data/sda/tzy/NA12878/GRCh38_reference_genome
wget -c https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa.fai /home/data/sda/tzy/NA12878/GRCh38_reference_genome
wget -c https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa.pac /home/data/sda/tzy/NA12878/GRCh38_reference_genome
wget -c https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa.sa /home/data/sda/tzy/NA12878/GRCh38_reference_genome
wget -c https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/README.20150309.GRCh38_full_analysis_set_plus_decoy_hla /home/data/sda/tzy/NA12878/GRCh38_reference_genome
- 使用samtools检测文件完整性,如果文件完整则不输出内容:
> samtools quickcheck NA12878.bam
NA12878.bam was missing EOF block when one should be present.
2. 转化
转成bam文件
#!/bin/bash
samtools view -@ 120 -T /home/data/sda/tzy/NA12878/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa -b /home/data/sda/tzy/NA12878/NA12878.alt_bwamem_GRCh38DH.20150718.CEU.low_coverage.cram -o /home/data/sda/tzy/NA12878/NA12878.bam
3. 计算深度
#!/bin/bash
cd /home/data/sda/tzy/NA12878/
samtools depth -@ 120 -a /home/data/sda/tzy/NA12878/NA12878.bam > /home/data/sda/tzy/NA12878/depth.txt
awk '{sum += $3} END {print "Average depth: ", sum/NR}' /home/data/sda/tzy/NA12878/depth.txt > /home/data/sda/tzy/NA12878/ave_depth.txt
注意这个深度的计算应该是针对sorted的bam文件。
4. 生成bai文件
#!/bin/bash
samtools index -@ 120 /home/data/sda/tzy/NA12878/NA12878.bam
5. 根据工具需要处理
freec需要每个染色体的fasta文件
for i in `cut -f1 GRCh38_full_analysis_set_plus_decoy_hla.fa.fai`;do samtools faidx GRCh38_full_analysis_set_plus_decoy_hla.fa $i > ./chr/$i.fa;done
samtools faidx GRCh38_full_analysis_set_plus_decoy_hla.fa HLA-C\*02:11 > chr/HLA-C*02:11.fa
# 不知道什么情况,这个区域只能手动生成,总是识别不出来
6. 判断单端还是双端测序
返回0就是单端,1就是双端
samtools view -h NA12878.bam | head -n 100 | samtools view -c -f 1