mapthis.clean in genetic map construction tutorial

[James E Specht]

Karl:  My students have been working their way thru your nicely documented tutorial for Genetic map construction with R/qtl, using your supplied command code for this tutorial.  However, one student who followed up on the mapthis.clean set of cmds (see below) with his own extra follow-up plot.geno command indicated to me that the eight genotypes did not seem to have been cleaned up (i.e., converted to NA or missing) in the mapthis.clean object.  I did note to the student that the solid or open symbols red-boxed in the 1st plot did disappear in the 2nd plot (in some cases replaced by an x symbol), but red boxes were still present in the 2nd plot at the same positions.  Can you offer us some insight or clarification on this?  Thanks!  - Jim  


> ### chunk number 122: plotgeno eval=FALSE
> ###################################################
> ## #line 1920 "geneticmaps.Rnw"
> ## plot.geno(mapthis, chr=1, ind=toperr$id[toperr$chr==1],
> ##           cutoff=6, include.xo=FALSE)
>
>
> ###################################################
> ### chunk number 123: plotgenoplot
> ###################################################
> #line 1927 "geneticmaps.Rnw"
> par(mar=c(4.1,4.1,0.6,0.6), las=1, cex.axis=0.9)
> plot.geno(mapthis, chr=1, ind=toperr$id[toperr$chr==1], main="", cex=0.8,
+           include.xo=FALSE, cutoff=6)
>
>
> ###################################################
> ### chunk number 124: dropgenotypes
> ###################################################
> #line 1953 "geneticmaps.Rnw"
> mapthis.clean <- mapthis
> for(i in 1:nrow(toperr)) {
+   chr <- toperr$chr[i]
+   id <- toperr$id[i]
+   mar <- toperr$marker[i]
+   mapthis.clean$geno[[chr]]$data[mapthis$pheno$id==id, mar] <- NA
+ }
>
>
> plot.geno(mapthis.clean, chr=1, ind=c("id200","id36","id236","id217","id235","id261","id87","id72"),cutoff=6,include.xo=FALSE)

[Karl Broman]

You need to first re-run calc.errorlod. plot.geno(mapthis.clean, ...) is looking at the previously calculated error LOD scores.

> mapthis.clean <- calc.errorlod(mapthis.clean, error.prob=0.005)
> top.errorlod(mapthis.clean, cutoff=6)
No errorlods above cutoff.

karl

> --
> You received this message because you are subscribed to the Google Groups "R/qtl discussion" group.
> To post to this group, send email to rqtl...@googlegroups.com.
> To unsubscribe from this group, send email to rqtl-disc+...@googlegroups.com.
> For more options, visit this group at http://groups.google.com/group/rqtl-disc?hl=en.

[e.b...@tum.de]

Dear Karl,

I have a similar problem. When I use chunk #124: dropgenotypes as shown above with my data (cross object gtypes_bin4_1to10), it seems that converting the erroneous genotypes to NA does not work.

errlod <- calc.errorlod(gtypes_bin4_1to10, error.prob=0.01, map.function=c("kosambi"), version=c("new"))
toperr <- top.errorlod(errlod, cutoff=3, msg=TRUE)
dim(toperr)
# [1] 3575    4

head(toperr)
# chr  id      marker errorlod
# 1  L4 154 M123 6.887229
# 2  L7  36 M234 5.285263
# only 2 rows shown here

# clean up the data:
gtypes_bin4_1to10.clean <- gtypes_bin4_1to10
  for(i in 1:nrow(toperr)) {
    chr <- toperr$chr[i]
    id <- toperr$id[i]
    mar <- toperr$marker[i]
    gtypes_bin4_1to10.clean$geno[[chr]]$data[gtypes_bin4_1to10$pheno$id==id, mar] <- NA
  }

I checked it by rerunning calc.errorlod on the *.clean object, but it still gives the same list of genotypes with top.errorlod (total of 3575 rows in the object).

errlod2 <- calc.errorlod(gtypes_bin4_1to10.clean, error.prob=0.01, map.function=c("kosambi"), version=c("new"))
toperr2 <- top.errorlod(errlod2, cutoff=3, msg=TRUE)
dim(toperr2)
# [1] 3575    4

My original data of 182 F2 plants and 8869 SNPs does not have any missing data. If I compare with the *.clean object, the data are exactly the same - no NAs.

Any suggestions what I could check or change are very welcome.

Best regards
Eva

 p.s.: I am using qtl version 1.41-6, R-3.4.2 and RStudio v 1.1.383.

[Karl Broman]

I can't tell what's going wrong. I'd break it down into smaller pieces. For example, look at the very first item in the loop: are chr, id, and mar taking the correct values?

Personally, I no longer try to remove genotypes that appear to be in error, but rather just use error.prob=0.01 or =0.002 in calc.genoprob() and est.map().

karl

> To unsubscribe from this group and stop receiving emails from it, send an email to rqtl-disc+...@googlegroups.com.

> To post to this group, send email to rqtl...@googlegroups.com.

> Visit this group at https://groups.google.com/group/rqtl-disc.
> For more options, visit https://groups.google.com/d/optout.

[e.b...@tum.de]

Dear Karl,

the problem was, that in the id column of toperr there were consecutive numbers for the genotypes, but in the cross object I had the "real" genotype IDs (like line1, line2 etc.). I changed the genotype IDs to consecutive numbers, too, then it worked.

Best regards
Eva