error relate to R

Question

[rheab1230]

I am running a R script to update vcf file with rsid of snp file.I am getting error.

Error in `$<-.data.frame`(`*tmp*`, variant_id, value = "snp__") :
  replacement has 1 row, data has 0
Calls: $<- -> $<-.data.frame
Execution halted

Can anyone please help me with this?

[ATpoint]

Hi, this is not reproducible. Both code and data are missing. Please provide it with a minimum working example.

[rheab1230]

## this is the R script:
## called from the matching .sh script (which passes in the parameters)


# getting file name as passed from the bash script
passed_args<-commandArgs(trailingOnly = TRUE)

# limits the snp reference file to only the chrom of interest
chr_num<-passed_args[3]

# variable for the snp151 reference file - expecting at least a CHROM, POS, and rsid columns in the snp ref file
snp_file<-passed_args[4]

# imports the reference data file for the SNPs (old vcf "ID" column contains snp_22_30192384; new "ID" should be the rsid)
snp_ref<-read.delim(snp_file, header=TRUE, sep=" ")
snp_ref<-snp_ref[which(snp_ref$CHROM==chr_num),]
snp_ref$variant_id<-paste("snp_",snp_ref$CHROM, "_",  snp_ref$POS, sep="")
snp_ref$CHROM<-as.integer(snp_ref$CHROM)
paste("Chr ", chr_num, "SNPs found: ", nrow(snp_ref), sep="")

#get from args[] the .vcf file
vcf_file_name<- passed_args[1]

#find the first line of actual data
header_rows_skip<-as.integer(passed_args[2])-1

#imports the file into a table without the metadata lines
body_only_VCF<-read.delim(vcf_file_name, header = TRUE, skip=header_rows_skip, sep="\t")
colnames(body_only_VCF)[1]<-"CHROM"
paste("orig VCF file rows ", nrow(body_only_VCF))

#removes anything where the reference or alternate allele is not G/C/T/A (i.e. narrows down vcf working file to just SNPs)
body_only_VCF<-body_only_VCF[which(body_only_VCF$REF %in% c("A","T","C","G") & body_only_VCF$ALT %in% c("A","T","C","G")),]
paste("VCF filtered to only SNPS ", nrow(body_only_VCF))

#merges the VCF file with records that have a matching lookup in the snp_ref table
final_VCF<-merge(body_only_VCF,snp_ref, all = FALSE )
paste("Merge complete - matched VCF records by CHROM and POS", nrow(final_VCF))

#get final file column order
VCF_col_names<-names(body_only_VCF)

#trying to free up some memory
rm(body_only_VCF)

#udpated the ID column and removes the columns from snp_ref
final_VCF$ID<-final_VCF$rsid
final_VCF<-final_VCF[,VCF_col_names]

#rename first col to #CHROM to match vcf file formatting
colnames(final_VCF)[1]<-"#CHROM"

#write final file to the vcf_directory
final_file=paste(substr(vcf_file_name, 1, nchar(vcf_file_name)-4), "_updated2.vcf", sep="")
write.table(final_VCF, file = final_file, sep="\t", row.names = FALSE, quote=FALSE)

[rheab1230]

What the script does is it takes rsid from snp file and put that in ID coloumn in vcf file.

[Kevin Blighe]

Thanks for showing the code; however, this doesn't help due to the fact that we cannot have certainty which line produces the error. You have neither provided any test input data; so, we cannot run the code in its entirety.

Please tell us which line produces the error. Also show the output of str() run on the objects that are used in the command that produces the error.

[rheab1230]

This is the line that is showing the error: snp_ref$variant_id<-paste("snp_",snp_ref$CHROM, "_", snp_ref$POS, sep="")

[rheab1230]

I am not able to attach the input data file. Can you please tell me how to attach the input file?

[Kevin Blighe]

I see, what is the output of str(snp_ref) just before you run that line of code?

[rheab1230]

str(snp_ref)
spec_tbl_df [9 x 5] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ CHROM   : num [1:9] 1 1 1 1 1 1 1 1 1
 $ POS     : num [1:9] 10038 10042 10050 10054 10062 ...
 $ chromEnd: num [1:9] 10039 10043 10051 10055 10063 ...
 $ rsid    : chr [1:9] "rs978760828" "rs1008829651" "rs1052373574" "rs892501864" ...
 $ X5      : chr [1:9] NA NA NA NA ...
 - attr(*, "problems")= tibble [9 x 5] (S3: tbl_df/tbl/data.frame)
  ..$ row     : int [1:9] 1 2 3 4 5 6 7 8 9
  ..$ col     : chr [1:9] NA NA NA NA ...
  ..$ expected: chr [1:9] "5 columns" "5 columns" "5 columns" "5 columns" ...
  ..$ actual  : chr [1:9] "4 columns" "4 columns" "4 columns" "4 columns" ...
  ..$ file    : chr [1:9] "'C:/Users/Desktop/variant_mapping_small.txt'" "'C:/Users/Desktop/variant_mapping_small.txt'" "'C:/Users/Desktop/variant_mapping_small.txt'" "'C:/Users/Desktop/variant_mapping_small.txt'" ...
 - attr(*, "spec")=
  .. cols(
  ..   CHROM = col_double(),
  ..   POS = col_double(),
  ..   chromEnd = col_double(),
  ..   rsid = col_character(),
  ..   X5 = col_character()
  .. )

[Kevin Blighe]

Oh, your object is a weird format.. spec_tbl_df ... please change it to a data frame via data.frame(snp_ref)

[rheab1230]

Okay, I will do that

[rheab1230]

I have one more question. When I am mapping vcf file and dbsnp file so that vcf file can contain rsid, my vcf file size is getting reduced.

[Kevin Blighe]

That makes sense. Please view the VCF via bcftools view my.vcf | more, and then you'll observe why (not all of your variants will have an rs ID).

[rheab1230]

The dbSNP file is a file coming from the UCSC Genome Browser that would have a list of all the known SNPs with their rsid. The GEUVADIS .vcf files I used were older and did not have an appropriate rsid in the ID column of that file. So I made a script to update .vcf file to assign them an appropriate modern rsid. However, the record count had dropped by a lot from the original .vcf file. Basically SNPs are getting lost from the .vcf file when merged with the dbSNP file.

[Kevin Blighe]

It could be variants in GEUVADIS that have no known dbSNP rs ID, I suppose