R/QTL for biparental population

[Pummi Singh]

Hi,

I am trying to generate a genetic map for a biparental fungal population (haploid). I am almost brand new to coding in R. I am stuck at the following code for recombination frequencies. Any help is appreciated. Thanks.

I can share the code before the following code if required. 

rf.(mapthis) <- pull.rf(mapthis)

lod.(mapthis) <- pull.rf(mapthis, what="lod")

plot(as.numeric(rf), as.numeric(lod), xlab="Recombination fraction", ylab="LOD score")

Error in `.rowNamesDF<-`(x, value = value) : duplicate 'row.names' are not allowed

8. stop("duplicate 'row.names' are not allowed")

7. `.rowNamesDF<-`(x, value = value)

6. `row.names<-.data.frame`(`*tmp*`, value = value)

5. `row.names<-`(`*tmp*`, value = value)

4. `rownames<-`(`*tmp*`, value = unlist(lapply(map, names)))

3. map2table(pull.map(cross, as.table = FALSE))

2. pull.map(cross, as.table = TRUE)

1. pull.rf(mapthis)

[Karl Broman]

It looks like you may have duplicate marker names.

Type 'summary(mapthis)' and you should get a warning about it.

karl

--
You received this message because you are subscribed to the Google Groups "R/qtl discussion" group.
To unsubscribe from this group and stop receiving emails from it, send an email to rqtl-disc+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/rqtl-disc/6e205d5d-f6cd-41ec-bd94-b6d3ef11ce40n%40googlegroups.com.

[Pummi Singh]

Hi,

Yes, I used drop duplicates and now the code is working.

Thank you Dr Broman.

Pummi

You received this message because you are subscribed to a topic in the Google Groups "R/qtl discussion" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/rqtl-disc/rjM-1NGGze0/unsubscribe.
To unsubscribe from this group and all its topics, send an email to rqtl-disc+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/rqtl-disc/CAEwE18H10vSrmDJtfUBa9eEGPHhYxiTDMAWR4aGMZq%3DQeC2ZbQ%40mail.gmail.com.

[Pummi Singh]

Hi,

I have mapped genotype frequencies by individual. Since I am working with a haploid organism, I only have AA and BB genotypes. One individual has genotype identical to BB (100%) and I want to drop this individual. Is it possible to drop an individual itself in R/QTL or do I drop markers for this individual? What function should I use to do so?

Thanks

[Karl Broman]

Use the function subset(). Something like:

  mycross_subset <- subset(mycross, ind = (pB < 1))

There's an example like this at the bottom of page 100 in the R/qtl book (https://rqtl.org/book), in chapter 4.

It's available as a sample chapter: https://rqtl.org/book/rqtlbook_ch04.pdf; see page 27 of the PDF.

karl

[Pummi Singh]

Thank you Dr.Broman.

To view this discussion on the web visit https://groups.google.com/d/msgid/rqtl-disc/12eaf5a8-9f81-40ff-b3be-18d790153123n%40googlegroups.com.

[Pummi Singh]

Hi,

I am trying to run ordermarkers on a linkage group with 438 markers. It keeps giving me the error "cannot allocate vector of size xx Gb", where xx depends on what system I am running it on. Is this a memory issue? Is it too much for my system to handle 438 markers on a single linkage group or in general 438 markers is too much? 

Thanks

Pummi 

[Karl Broman]

Yes that would indicate a memory problem, but I can’t see hoe it would come up with only 438 markers. Can you provide more drtails: the exact code you’re using, and where it is that the error occurs?

karl

To view this discussion on the web visit https://groups.google.com/d/msgid/rqtl-disc/9796e01f-41b5-4f5c-be3d-77f8d9dcfd58n%40googlegroups.com.

[Pummi Singh]

Hi Dr.Broman,

I am trying to run the following after running all the previous steps:

I was able to run chr 2 with 290 markers.

```{r}
map <- orderMarkers(map, chr=1)
pull.map(map, chr=1)
```

The error is:

Error: cannot allocate vector of size 5.9 Gb

5.	rbind(orders, cbind(left, temp + i - 1, right))

4.	ripple(cross, chr = 1, window = window, verbose = FALSE)

3.	summary(ripple(cross, chr = 1, window = window, verbose = FALSE))

2.	orderMarkers.sub(cross, i, window = window, use.ripple = use.ripple, verbose = verbose.sub)

1.	orderMarkers(map, chr = 1)

To view this discussion on the web visit https://groups.google.com/d/msgid/rqtl-disc/7A7605D7-9D8C-4C82-827F-6A8DB1442018%40gmail.com.

[Karl Broman]

It's running out of memory in ripple(). I would either try use.ripple=FALSE or reduce the window size.

   map <- orderMarkers(map, chr=1, use.ripple=FALSE)

karl

[Pummi Singh]

Hi Dr.Broman,

I had a quick question about ripple function. Please find the document attached. In the first screenshot, certain orders do not converge (shown as "Didn't converge!" in red). I was wondering if that is a problem that needs attention.

In the second screenshot, different orders of the markers with LOD scores are present. It is my understanding that Order 1 is preferred over the Initial (also the rest of them) because of high LOD and lower chromosome length. Also, it appears that every consecutive order of markers is presented such that LODs are in decreasing order (because I am not able to see what data is in those 12 additional rows). Please let me know if this is true.

Thanks

Pummi 

[Karl Broman]

If you use the likelihood method in ripple(), it will re-estimate inter-marker distances with each marker order, using an iterative algorithm. The warning message is issued if that algorithm doesn't converge. You can increase the maxit argument to try to achieve convergence, but it probably doesn't matter and can be ignored, as the lack of convergence is mostly happening for marker orders with lower support.

The LOD scores are log10 likelihood relative to the initial order, so yes the one with the largest value is preferred. See the help file for ripple(): "the last two columns are LOD scores (log base 10 likelihood ratios) comparing each order to the initial order and the estimated chromosome length for the given order. Positive LOD scores indicate that the alternate order has more support than the original."

karl

[Pummi Singh]

Great. Thank you Dr.Broman.

To view this discussion on the web visit https://groups.google.com/d/msgid/rqtl-disc/6544eac6-aba7-4126-8d65-ecc45ec3495en%40googlegroups.com.

[Pummi Singh]

Hi Dr.Broman,

What do you suggest in situations when LOD score is high after running ripple but the chromosome length also increases? I am particularly concerned when chromosome length increases by more than 1 centiMorgan. Please find the screenshot of one of my example codes' attached.

Thank you. 

[Karl Broman]

Without any other information about physical order of the markers, I'd go with the order with highest likelihood, even if the estimated chromosome length increases slightly.

karl

[Pummi Singh]

Thank you Dr.Broman.

I am getting the following error while running the code for estimating genotyping error rate. Please let me know how to fix this.

Thanks

 

[Pummi Singh]

I am attaching a word doc with the code and error because I am not sure if what I pasted is visible or not. 

Thanks

[Karl Broman]

I would pull it apart and try to figure out what's going wrong.

Look at the output of sapply(tempmap, attr, "loglik") which should be a vector of log likelihood values.

karl

[Pummi Singh]

Hi Dr.Broman,

Attached are the results that I got. I have 17 linkage groups of which groups 15 to 17 have a single marker each. I apologize for several questions, I am very new to R. 

Thanks

[Karl Broman]

Drop the groups with just one marker. They don't have a log likelihood and so are messing things up.

karl

To view this discussion on the web visit https://groups.google.com/d/msgid/rqtl-disc/df2430a5-f0e3-4fa3-8337-669af1564ad6n%40googlegroups.com.

[Pummi Singh]

Does that mean dropping the markers associated with these groups?

Thanks

[Karl Broman]

I was thinking that you'd drop the chromosomes with subset(), or get around those chromosomes when you try to pull the log likelihood values.

Since the problem is that sapply( ) is returning a list, due to some chromosomes returning NULL, you could change that line to have unlist(lapply(  )), as that will strip of the NULL values:

    sum(unlist(lapply(tempmap, attr, "loglik"))) 

karl

PS It would be a bit easier to respond/work with your code if you copy-pasted it into your message rather than adding a screen shot image.

[Pummi Singh]

Thanks Dr.Broman. Sorry about the screenshots. I will make sure that I have the code and error pasted instead.

To view this discussion on the web visit https://groups.google.com/d/msgid/rqtl-disc/8e54eefb-6882-4f4e-8956-47486ab9525en%40googlegroups.com.

[Pummi Singh]

Hi Dr.Broman,

I deleted a marker thinking that it might improve my results but that did not happen. I want to insert that marker back on its linkage group. I have so far not been successful at finding a function to insert a deleted marker into a linkage group. Can you please suggest something.

Thank you.

Pummi

[Karl Broman]

There's not a way to put it back if you've deleted it from the cross. Rather than delete markers, I would move them to some "unlinked" linkage group or chromosome. To get the data back, you'll need to go back to the original data set.

karl

[Pummi Singh]

Ahh ok.. Thanks

To view this discussion on the web visit https://groups.google.com/d/msgid/rqtl-disc/ae8a6fb5-8b08-4137-be03-1afc95addf3an%40googlegroups.com.

[Pummi Singh]

Hi Dr.Broman,

I had to move certain markers from one linkage group to another and as a result, the first marker of the initial linkage group starts from a number other than 0 cM. I am not sure if this is fine. Please let me know if there is some way to fix this. An example is attached. I had to move the markers highlighted in orange to another linkage group. After moving the markers, I ran the following for the linkage group to and from which markers were moved. 

```{r}

genmap <- orderMarkers(genmap, chr=2)

pull.map(genmap, chr=2)

```

```{r}

rip2 <- ripple(genmap, chr=2, window=7)

summary(rip2)

```

```{r}

rip2lik <- ripple(genmap, chr=2, window=4, method="likelihood", error.prob=0.005)

summary(rip2lik)

```

I also used switch.order after ripple based on LOD values.

Please let me know.

Thanks

Pummi

[Karl Broman]

There's a function shiftmap() that can be used to shift all chromosomes so they start at 0.

karl

[Pummi Singh]

Thank you.

To view this discussion on the web visit https://groups.google.com/d/msgid/rqtl-disc/f83fab71-dba8-4d63-95ee-a2a74660a9c7n%40googlegroups.com.

[Pummi Singh]

Hi Dr.Broman,

How can I arrange my linkage groups in descending order of length? I had to make some new linkage groups and want to make sure that linkage group 1 is the longest followed by 2 and so on. 

Thank you

Pummi

[Karl Broman]

First pull out the map and get the lengths of the linkage groups. I'll use the cross 'fake.bc' as an example.

  library(qtl)

  data(fake.bc)

  map <- pull.map(fake.bc)

  L <- sapply(map, function(chr) diff(range(chr)))

Then reorder the chromosomes by reordering the components of the "geno" part of the cross.

  fake.bc$geno <- fake.bc$geno[order(L, decreasing=TRUE)]

If you want to change the names of the groups, you can change the names of the "geno" part of the cross

  names(fake.bc$geno) <- 1:nchr(fake.bc)

karl

[Pummi Singh]

Thanks so much Dr.Broman. 

To view this discussion on the web visit https://groups.google.com/d/msgid/rqtl-disc/e85b2288-5f5b-4a8d-99a1-d7892f71a85cn%40googlegroups.com.