Search
Question: DESeq2 results export padj sorted
0
gravatar for rbronste
4 months ago by
rbronste60
rbronste60 wrote:

When I do the following for some reason the output only includes ranges from a single chromosome or several, while if I look at the padj sorted GRanges in R it includes many chromosomes at the top of the padj list. Is there something I am doing wrong in the sorting and exporting of results? Thanks.

bn_group_male<-results(dds, lfcThreshold=1, altHypothesis="greater", format = c("GRanges"), independentFiltering=TRUE, contrast=c("group", "maleBB", "malevehicle"))

bn_group_male_padj.1 <- subset(bn_group_male, padj<.1)

bn_group_male_padj.1_sorted = bn_group_male_padj.1[order(bn_group_male_padj.1$padj<.1), ]

write.table(head(bn_group_male_padj.1_sorted,1000),"bn_group_male_padj.1.bed", quote=F, sep="\t", row.names=F, col.names=T)
ADD COMMENTlink modified 4 months ago by Michael Love17k • written 4 months ago by rbronste60

You have: 

bn_group_male_padj.1_sorted = bn_group_male_padj.1[order(bn_group_male_padj.1$padj<.1), ]

I guess this is not what you meant. If you wanted to use "order", there shouldn't be "<.1".

ADD REPLYlink written 4 months ago by Simon Anders3.5k

Well basically in the end I want a BED file of the top 1000 ranges padj<.1, in this list there are many thousands however I just want to export the top thousand by padj.

ADD REPLYlink written 4 months ago by rbronste60

It doesn't work that they. (You should have indicated that you are a novice R user.)

You first subset, and then you sort:

bn_group_male_padj.1 <- subset(bn_group_male, padj<.1)
bn_group_male_padj.1_sorted = bn_group_male_padj.1[order(bn_group_male_padj.1$padj), ]
ADD REPLYlink written 4 months ago by Simon Anders3.5k

I am a fairly novice user however if you look above at the initial example I did in fact first subset and then sort. Thanks.

ADD REPLYlink written 4 months ago by rbronste60
1

I know. But you put the "<.1" into the parantheses after "order". So what happens is that the argument of order gets transformed into a vector of only TRUE and FALSE. The actual values of padj vanish. So, the ordering only ensures that FALSE comes before TRUE. 

Now, as you have already filtered, there is only TRUE, and the line with "order" does nothing. Everything stays in the order it was in before. And that was presumably ordered by chromosome. That's why the first 1000 items were all on chromosome 1.

What I just did is copy your two lines and remove the "<.1". 

ADD REPLYlink written 4 months ago by Simon Anders3.5k

Thanks for the help!

ADD REPLYlink written 4 months ago by rbronste60
1
gravatar for Michael Love
4 months ago by
Michael Love17k
United States
Michael Love17k wrote:

See rtracklayer for writing GRanges to BED file.

ADD COMMENTlink written 4 months ago by Michael Love17k

Thanks for the suggestion!

While I do get the bed file I am after, it seems to be sorted by chromosome as it only gives me chr1 ranges and not other. Not sure why since I am doing the following:

bn_group_male_padj.1_sorted = bn_group_male_padj.1[order(bn_group_male_padj.1$padj<.1), ]

export.bed(head(bn_group_male_padj.1_sorted,1000),"bn_group_male_padj.1.bed")​

So it seems like I am sorting by padj and all I basically want here is ranges on any chromosome that fall into the top 1000 by padj. Not exactly sure what I am doing wrong, thanks again.

ADD REPLYlink modified 4 months ago • written 4 months ago by rbronste60

I can't tell when you have a data.frame and when you have a GRanges. Does subset() give back a GRanges? I'd be careful about when you've got which class of object. You can definitely use [] and which(res$padj < .1) to operate on a GRanges. And [] with order(res$padj) should work as well.

ADD REPLYlink written 4 months ago by Michael Love17k

Yes they are GRanges so not sure why its not ordering them properly. 

ADD REPLYlink written 4 months ago by rbronste60
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.2.0
Traffic: 135 users visited in the last hour