Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • Bedtools intersectBed with vcf format

    Since I could not find any documentation on how Bedtools intersectBed works with vcf format, I'm posting my question here:

    I'm intersecting a vcf file from mpileup and bcftools with dbSNP132 that is in vcf format:

    intersectBed -a mutations_chr8.vcf -b dbSNP132_chr8.vcf -wao > mutations_chr8_dbSNP132.txt

    In my example below I expected to see only dbSNP entries that have the exactly same coordinates, ie.e chr8 112454855, but I also get hits for chr8 112454865 and chr8 112454876:

    Code:
    chr8	112454855	.	Tacacacacacacacacacacaca	TACacacacacacacacacacacaca	35.5	.	INDEL;DP=24;AF1=0.5;CI95=0.5,0.5;DP4=0,3,2,2;MQ=46;PV4=0.43,0.3,0.34,1	PL:GT:GQ	73,0,78:0/1:75	chr8	112454855	rs76365113	T	A	.	.	dbSNPBuildID=131;VP=050000000001070008000100;WGT=1;VC=SNP;VLD;G5A;G5;KGPilot1	1
    chr8	112454855	.	Tacacacacacacacacacacaca	TACacacacacacacacacacacaca	35.5	.	INDEL;DP=24;AF1=0.5;CI95=0.5,0.5;DP4=0,3,2,2;MQ=46;PV4=0.43,0.3,0.34,1	PL:GT:GQ	73,0,78:0/1:75	chr8	112454855	rs67322031	A	AAC	.	.	dbSNPBuildID=130;VP=050000000001000000000200;WGT=1;VC=INDEL	1
    chr8	112454855	.	Tacacacacacacacacacacaca	TACacacacacacacacacacacaca	35.5	.	INDEL;DP=24;AF1=0.5;CI95=0.5,0.5;DP4=0,3,2,2;MQ=46;PV4=0.43,0.3,0.34,1	PL:GT:GQ	73,0,78:0/1:75	chr8	112454855	rs71851070	TAC	T	.	.	dbSNPBuildID=130;VP=050000000001000000000200;WGT=1;VC=INDEL	3
    chr8	112454855	.	Tacacacacacacacacacacaca	TACacacacacacacacacacacaca	35.5	.	INDEL;DP=24;AF1=0.5;CI95=0.5,0.5;DP4=0,3,2,2;MQ=46;PV4=0.43,0.3,0.34,1	PL:GT:GQ	73,0,78:0/1:75	chr8	112454865	rs34424884	CAC	C	.	.	dbSNPBuildID=126;VP=050000000001000000000200;WGT=1;VC=INDEL	3
    chr8	112454855	.	Tacacacacacacacacacacaca	TACacacacacacacacacacacaca	35.5	.	INDEL;DP=24;AF1=0.5;CI95=0.5,0.5;DP4=0,3,2,2;MQ=46;PV4=0.43,0.3,0.34,1	PL:GT:GQ	73,0,78:0/1:75	chr8	112454876	rs36066513	ACA	A	.	.	dbSNPBuildID=126;VP=050100000001030100000200;WGT=1;VC=INDEL;SLO;G5A;G5;GNO	3
    This is correct and useful considering that indel coordinates are often a bit shifted. However, I'd like to how Bedtools determines the interval in which to look for the intersection. I guess it determines the length of the reference allele the -a file, but does it also do the same for the -b file? And how could I restrict the intersection to the exact position instead of the interval?

    Thanks in advance for any help - and for Bedtools in general

    Barbara

  • #2
    Hi Barbara

    If you generate a bed file from your VCF file you can define intervals for one base like this:

    Code:
    chr8	112454855 112454856 Tacacacacacacacacacacaca	TACacacacacacacacacacacaca	35.5	.	INDEL;DP=24;AF1=0.5;CI95=0.5,0.5;DP4=0,3,2,2;MQ=46;PV4=0.43,0.3,0.34,1	PL:GT:GQ	73,0,78:0/1:75
    Using shell this can be done by:

    Code:
    cat mutations_chr8.vcf | awk '{print $1"\t"$2"\t"$2+1"\t"$0}' | cut -f -3,6- > mutations_chr8.bed
    
    cat dbSNP132_chr8.vcf  | awk '{print $1"\t"$2"\t"$2+1"\t"$0}' | cut -f -3,6- > dbSNP132_chr8.bed
    Then using intersectBed will match exact position.

    Comment


    • #3
      Good trick. Just a minor correction, the VCF pos is 1-based and BED pos is 0-based, so the awk part should be:
      awk '{print $1"\t"$2-1"\t"$2"\t"$0}'

      Comment

      Latest Articles

      Collapse

      • seqadmin
        Exploring the Dynamics of the Tumor Microenvironment
        by seqadmin




        The complexity of cancer is clearly demonstrated in the diverse ecosystem of the tumor microenvironment (TME). The TME is made up of numerous cell types and its development begins with the changes that happen during oncogenesis. “Genomic mutations, copy number changes, epigenetic alterations, and alternative gene expression occur to varying degrees within the affected tumor cells,” explained Andrea O’Hara, Ph.D., Strategic Technical Specialist at Azenta. “As...
        07-08-2024, 03:19 PM
      • seqadmin
        Exploring Human Diversity Through Large-Scale Omics
        by seqadmin


        In 2003, researchers from the Human Genome Project (HGP) announced the most comprehensive genome to date1. Although the genome wasn’t fully completed until nearly 20 years later2, numerous large-scale projects, such as the International HapMap Project and 1000 Genomes Project, continued the HGP's work, capturing extensive variation and genomic diversity within humans. Recently, newer initiatives have significantly increased in scale and expanded beyond genomics, offering a more detailed...
        06-25-2024, 06:43 AM

      ad_right_rmr

      Collapse

      News

      Collapse

      Topics Statistics Last Post
      Started by seqadmin, Yesterday, 07:20 AM
      0 responses
      24 views
      0 likes
      Last Post seqadmin  
      Started by seqadmin, 07-16-2024, 05:49 AM
      0 responses
      38 views
      0 likes
      Last Post seqadmin  
      Started by seqadmin, 07-15-2024, 06:53 AM
      0 responses
      44 views
      0 likes
      Last Post seqadmin  
      Started by seqadmin, 07-10-2024, 07:30 AM
      0 responses
      41 views
      0 likes
      Last Post seqadmin  
      Working...
      X