Skip to content

Instantly share code, notes, and snippets.

@fo40225
Last active October 21, 2021 04:31
Show Gist options
  • Save fo40225/8cfd1a6d84d7a1e641421ed1bf01c624 to your computer and use it in GitHub Desktop.
Save fo40225/8cfd1a6d84d7a1e641421ed1bf01c624 to your computer and use it in GitHub Desktop.
create clinical whole exome sequencing bed file
wget http://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_15/gencode.v15.annotation.gtf.gz
wget http://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_21/gencode.v21.chr_patch_hapl_scaff.annotation.gtf.gz
gunzip gencode.v15.annotation.gtf.gz
gunzip gencode.v21.chr_patch_hapl_scaff.annotation.gtf.gz
grep exon gencode.v15.annotation.gtf \
|awk '{printf("%s\t%s\t%s\n",$1,int($4)-1,$5);}' |sort -T . -k1,1 -k2,2n | bedtools merge > gencode.37.bed
grep exon gencode.v21.chr_patch_hapl_scaff.annotation.gtf \
|awk '{printf("%s\t%s\t%s\n",$1,int($4)-1,$5);}' |sort -T . -k1,1 -k2,2n | bedtools merge > gencode.38.bed
# edit bed, remove unnecessary chromosome
bedtools sort -faidx bundle/hg19/ucsc.hg19.fasta.fai -i gencode.37.bed >WES_GENCODE.37.bed
bedtools sort -faidx bundle/hg38/Homo_sapiens_assembly38.fasta.fai -i gencode.38.bed >WES_GENCODE.38.bed
sed -r -i 's/^chr//g' WES_GENCODE.37.bed
sed -r -i 's/^chr//g' WES_GENCODE.38.bed
wget --no-passive-ftp -O clinvar_20210626.37.vcf.gz ftp://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh37/archive_2.0/2021/clinvar_20210626.vcf.gz
wget --no-passive-ftp -O clinvar_20210626.38.vcf.gz ftp://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/archive_2.0/2021/clinvar_20210626.vcf.gz
gunzip clinvar_20210626.37.vcf.gz
gunzip clinvar_20210626.38.vcf.gz
bcftools query -i 'CLNSIG == "Pathogenic" || CLNSIG == "Likely_pathogenic"' -o 37.clinvar.bed \
-f '%CHROM\t%POS0\t%POS\n' clinvar_20210626.37.vcf
bcftools query -i 'CLNSIG == "Pathogenic" || CLNSIG == "Likely_pathogenic"' -o 38.clinvar.bed \
-f '%CHROM\t%POS0\t%POS\n' clinvar_20210626.38.vcf
bcftools query -i 'CLASS == "DM" || CLASS == "DM?"' -o 37.hgmd.bed \
-f '%CHROM\t%POS0\t%POS\n' HGMD_Pro_2021.2_hg19.vcf
bcftools query -i 'CLASS == "DM" || CLASS == "DM?"' -o 38.hgmd.bed \
-f '%CHROM\t%POS0\t%POS\n' HGMD_Pro_2021.2_hg38.vcf
cat WES_GENCODE.37.bed 37.clinvar.bed 37.hgmd.bed > 37.bed
cat WES_GENCODE.38.bed 38.clinvar.bed 38.hgmd.bed > 38.bed
sort -k1,1 -k2,2n 37.bed > 37.sorted.bed
sort -k1,1 -k2,2n 38.bed > 38.sorted.bed
bedtools merge -i 37.sorted.bed > 37.merged.bed
bedtools merge -i 38.sorted.bed > 38.merged.bed
# edit bed, remove unnecessary chromosome
cp 37.merged.bed 19.merged.bed
sed -r -i 's/^/chr/g' 19.merged.bed
sed -r -i 's/^/chr/g' 38.merged.bed
bedtools sort -faidx bundle/b37/human_g1k_v37_decoy.fasta.fai -i 37.merged.bed >wes.b37.bed
bedtools sort -faidx bundle/hg19/ucsc.hg19.fasta.fai -i 19.merged.bed >wes.hg19.bed
bedtools sort -faidx bundle/hg38/Homo_sapiens_assembly38.fasta.fai -i 38.merged.bed >wes.hg38.bed
gatk IndexFeatureFile -I wes.b37.bed
gatk IndexFeatureFile -I wes.hg19.bed
gatk IndexFeatureFile -I wes.hg38.bed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment