bedtools intersect error: Invalid record in file

Question

[Michal Nevo]

Hello to all

I am trying to run bedtools intersect with vcf file and a bed file (my goal is to add the depth data to my VCF) I get an error running this command:

bedtools intersect -a depth.bed -b fish.vcf -wa -wb > $out

The error:

Error: Invalid record in file depth.bed Record is 
NC_048323.1 1   0   0   0   0   1   0   00  0   0

I have a bed file I created with samtools depth

With the command : samtools depth -aa -H -f $list > $out_file

But you can see my bed file is fine:

NC_048323.1     1       0       0       0       0       1       0       0       0       0       0
NC_048323.1     2       0       0       0       1       1       0       0       0       0       0
NC_048323.1     3       0       0       0       1       1       0       0       0       0       0
NC_048323.1     4       0       0       0       2       2       0       0       0       0       0
.......
NC_022453.1     16788   17      30      12      25      20      12      18      13      21      17
NC_022453.1     16789   13      26      7       21      16      9       15      10      11      11
NC_022453.1     16790   5       13      5       11      10      6       6       8       8       6

Also, I checked if my bed file is tab-delimited:

NC_048323.1^I1^I0^I0^I0^I0^I1^I0^I0^I0^I0^I0$
NC_048323.1^I2^I0^I0^I0^I1^I1^I0^I0^I0^I0^I0$
NC_048323.1^I3^I0^I0^I0^I1^I1^I0^I0^I0^I0^I0$
NC_048323.1^I4^I0^I0^I0^I2^I2^I0^I0^I0^I0^I0$

I have the vcf file which I created with samtools depth

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  M1      M3      M5      M7      M9      F1.sorted.bamF2.sorted.bam    F4.sorted.bam   F6.sorted.bam   F8.sorted.bam
NC_048323.1     4135578 .       A       G       1473.64 .       DPB=25;EPPR=0;GTI=0;MQMR=0;NS=1;NUMALT=1;ODDS=18.7016;PAIREDR=0;PQR=0;PRO=0;QR=0;RO=0;RPPR=0;SRF=0;SRP=0;SRR=0;DP=186;AB=0;ABP=0;AF=1;AO=37;CIGAR=1X;DPRA=0;EPP=3.06899;LEN=1;MEANALT=1;MQM=60;PAIRED=0.972973;PAO=0;PQA=0;QA=1288;RPL=15;RPP=5.88603;RPR=22;RUN=1;SAF=16;SAP=4.47751;SAR=21;TYPE=snp;AN=20;AC=20   GT:DP:AD:RO:QR:AO:QA  ./././.:.:.:.:.:.:.     ./././.:.:.:.:.:.:.     ./././.:.:.:.:.:.:.     ./././.:.:.:.:.:.:.     ./././.:.:.:.:.:.:.   1/1/1/1:25:0,25:0:0:25:869      1/1/1/1:43:0,43:0:0:43:1498     1/1/1/1:48:0,48:0:0:48:1670     1/1/1/1:33:0,33:0:0:33:1167   1/1/1/1:37:0,37:0:0:37:1288
NC_048323.1     4465457 .       C       A       1244.59 .       DPB=27;EPPR=0;GTI=0;MQMR=0;NS=1;NUMALT=1;ODDS=19.8523;PAIREDR=0;PQR=0;PRO=0;QR=0;RO=0;RPPR=0;SRF=0;SRP=0;SRR=0;DP=163;AB=0;ABP=0;AF=1;AO=31;CIGAR=1X;DPRA=0;EPP=6.44263;LEN=1;MEANALT=1;MQM=60;PAIRED=0.935484;PAO=0;PQA=0;QA=1096;RPL=21;RPP=11.486;RPR=10;RUN=1;SAF=20;SAP=8.68415;SAR=11;TYPE=snp;AN=20;AC=20    GT:DP:AD:RO:QR:AO:QA  ./././.:.:.:.:.:.:.     ./././.:.:.:.:.:.:.     ./././.:.:.:.:.:.:.     ./././.:.:.:.:.:.:.     ./././.:.:.:.:.:.:.   1/1/1/1:27:0,27:0:0:27:942      1/1/1/1:30:0,30:0:0:30:1075     1/1/1/1:40:0,40:0:0:40:1414     1/1/1/1:35:0,35:0:0:35:1241   1/1/1/1:31:0,31:0:0:31:1096

What can I do?

Answer 1

[alex.zaccaron]

Your coverage file is not in bed format. The second and third columns of your $out_file must be start and end coordinates. Add a new second column that is equal to your current second column minus 1, with something like:

paste <(awk -v OFS="\t" '{print $1,$2-1}'  $out_file) <(cut -f 2-13 $out_file)

[Michal Nevo]

I just figured it out too! Thanks