通过VCF文件制作定制化的参考基因组

$ /gatk_software_path/gatk FastaAlternateReferenceMaker \
   -R /ref_genome_path/reference-genome.fasta 
   -O /ref_genome_path/psuedo-genome.fa \
   -L input.intervals \
   -V /variant_path/input.vcf \

   [--snp-mask mask.vcf]

$ perEditor ref_genome.fa snps_indels.vcf allele indivicual personalized_genome.fa

$ perEditor chr1.fa variations.vcf mother 1 chr1_customized.fa

$ ls chr*fa
chr1.fa  chr2.fa  chr3.fa  chr4.fa
$ perEditor_ra rearrangements.vcf chr_length.bed centromeres.bed individual allele
$ perEditor_ra rearrangements_t.vcf chr_length_t.bed centromeres_t.bed 1 father

$ ls *new*fa
chr1_new.fa  chr2_new.fa  chr3_new.fa  chr4_new.fa

$ conda config --add channels r
$ conda config --add channels bioconda
$ conda create -n g2gtools python=2 jupyter ipykernel
$ source activate g2gtools
$ conda install -c kbchoi g2gtools

source activate g2gtools

# 创建映射两个基因组之间的碱基的chain文件
g2gtools vcf2chain -f ${REF} -i ${VCF_INDELS} -s ${STRAIN} -o ${STRAIN}/REF-to-${STRAIN}.chain
# 把snp添补到ref genome上
g2gtools patch -i ${REF} -s ${STRAIN} -v ${VCF_SNPS} -o ${STRAIN}/${STRAIN}.patched.fa
# 将indels整合到snp添补后的基因组中
$ g2gtools transform -i ${STRAIN}/${STRAIN}.patched.fa -c ${STRAIN}/REF-to-${STRAIN}.chain -o ${STRAIN}/${STRAIN}.fa
# 针对新的定制基因组创建定制基因注释
$ g2gtools convert -c ${STRAIN}/REF-to-${STRAIN}.chain -i ${GTF} -f gtf -o ${STRAIN}/${STRAIN}.gtf
$ g2gtools gtf2db -i ${STRAIN}/${STRAIN}.gtf -o ${STRAIN}/${STRAIN}.gtf.db

git clone https://github.com/abyzovlab/vcf2diploid.git
cd vcf2diploid
make

java -Xmx10g -jar vcf2diploid.jar -id sample_id -chr file1.fa file2.fa ... [-vcf file1.vcf file2.vcf ...] > logfile.txt
OPTIONS:
id    - (required) the ID of individual whose genome is being constructed (e.g., NA12878). The tool recognizes by this ID in the VCF file 
chr   - (required) FASTA file(s) of reference sequence(s) 
vcf   - (required) VCF4.0 file(s) containing variants from parents and the individual 
Xmx   - max memory allocation for JAVA. In this example, 10GB was allocated.
logfile.txt - stores the standard output produce from the run

About: Create consensus sequence by applying VCF variants to a reference fasta
       file. By default, the program will apply all ALT variants. Using the
       --sample (and, optionally, --haplotype) option will apply genotype
       (or haplotype) calls from FORMAT/GT. The program ignores allelic depth
       information, such as INFO/AD or FORMAT/AD.
Usage:   bcftools consensus [OPTIONS] <file.vcf>
Options:
    -c, --chain <file>         write a chain file for liftover
    -e, --exclude <expr>       exclude sites for which the expression is true (see man page for details)
    -f, --fasta-ref <file>     reference sequence in fasta format
    -H, --haplotype <which>    choose which allele to use from the FORMAT/GT field, note
                               the codes are case-insensitive:
                                   1: first allele from GT
                                   2: second allele
                                   R: REF allele in het genotypes
                                   A: ALT allele
                                   LR,LA: longer allele and REF/ALT if equal length
                                   SR,SA: shorter allele and REF/ALT if equal length
    -i, --include <expr>       select sites for which the expression is true (see man page for details)
    -I, --iupac-codes          output variants in the form of IUPAC ambiguity codes
    -m, --mask <file>          replace regions with N
    -o, --output <file>        write output to a file [standard output]
    -s, --sample <name>        apply variants of the given sample
Examples:
   # Get the consensus for one region. The fasta header lines are then expected
   # in the form ">chr:from-to".
   samtools faidx ref.fa 8:11870-11890 | bcftools consensus in.vcf.gz > out.fa

$ bcftools consensus <file.vcf> \
  --fasta-ref <file> \
  --iupac-codes \
  --output <file> \
  --sample <name>
 
# Apply variants present in sample "NA001", output IUPAC codes for hets
$ bcftools consensus -i -s NA001 -f in.fa in.vcf.gz > out.fa

# Create consensus for one region. The fasta header lines are then expected
# in the form ">chr:from-to".
$ samtools faidx ref.fa 8:11870-11890 | bcftools consensus in.vcf.gz -o out.fa

$ samtools mpileup -d8000 -q 20 -Q 10 -uf  REFERENCE.fasta Your_File.bam | bcftools call -c | vcfutils.pl vcf2fq  > OUTPUT.fastq
# d, --max-depth
# -q, -min-MQ Minimum mapping quality for an alignment to be used
# -Q, --min-BQ Minimum base quality for a base to be considered
$ sed -i -e '/^+/,/^\@/{/^+/!{/^\@/!d}}; /^+/ d; s/@/>/g' OUTPUT.fastq