Question

ROH and Manhattan plot

0

Entering edit mode

4 months ago

Star:) ▴ 50

 --file genQC --homozyg --homozyg-group --homozyg-het 0
--homozyg-density 100 --homozyg-kb 1000 --homozyg-snp 45 --homozyg-window-snp  45  --homozyg-gap 1000 --homozyg-window-het 1 --homozyg-window-threshold 0.05 --homozyg-window-missing 0 --cow --out genQ

Hi all I want to draw a Manhattan plot for ROH results, but I am not sure which output (.hom, .overlap, .summary) I need to use. Can you guide me?

Thanks

ROH Manhattan plink • 731 views

ADD COMMENT • link updated 4 months ago by bk11 ★ 3.0k • written 4 months ago by Star:) ▴ 50

1

Entering edit mode

You can used .hom.overlap and basically plot on X-axis SNPs genomic coordinate, and the Y-axis showing the frequency (%) of overlapping ROH shared among individuals.

ADD REPLY • link 4 months ago by bk11 ★ 3.0k

0

Entering edit mode

Hi I am confused on how and which columns of file I can use to draw a manhattan plot for x axis information? I need to use BP1, BP2, CHR?

> head(hom_overlap)

  POOL FID    IID PHE CHR                  SNP1               SNP2     BP1

1   S1  kk 919  -9  1   ngs-102610 bovinehd09652 816847


       BP2        KB NSNP NSIM GRP

1 600580 1783.734  249 149   1

ADD REPLY • link 4 months ago by Star:) ▴ 50

1

Entering edit mode

CHR & BP1 need to be used for X-axis.

ADD REPLY • link 4 months ago by bk11 ★ 3.0k

0

Entering edit mode

Thanks for your help. I tried to write its code, but I am unsure about its accuracy. I used IID, which shows individuals. Is it true for frequency%? Could you please guide me again?

library(CMplot)

hom_overlap <- read.table("ks6.hom.overlap", header = TRUE)

roh_frequency <- aggregate(hom_overlap$IID, 
                           by = list(CHR = hom_overlap$CHR, BP1 = hom_overlap$BP1), 
                           FUN = function(x) length(unique(x)))

colnames(roh_frequency) <- c("CHR", "BP1", "individual_count")

total_individuals <- length(unique(hom_overlap$IID)) roh_frequency$frequency <- (roh_frequency$individual_count / total_individuals) * 100

hom_overlap_with_freq <- merge(hom_overlap, roh_frequency[, c("CHR", "BP1", "frequency")], by = c("CHR", "BP1"))

output_df <- hom_overlap_with_freq[, c("SNP1", "CHR", "BP1", "frequency")] colnames(output_df) <- c("SNP", "chr", "pos", "frequency")


write.table(output_df, "hom_overlap_with_frequency.txt", sep = "\t", row.names = FALSE, quote = FALSE)


CMplot(output_df, plot.type = "m", 
       band = 0.5, LOG10 = FALSE,  ylab = "Frequency", 
       threshold = NULL,  threshold.lty = 2, 
       threshold.lwd = 1,   threshold.col = "red", 
       amplify = TRUE,    width = 14, 
       height = 6,  signal.col = NULL,  chr.den.col = NULL, 
       file = "jpg", file.name = "manhattan_plot", 
       dpi = 300,   file.output = TRUE,    verbose = TRUE,   cex = 0.8)

Thank you very much.

ADD REPLY • link 4 months ago by Star:) ▴ 50

1

Entering edit mode

It looks good to me. I was able to generate ROH Manhattan plot using your the above code and my test data.

ADD REPLY • link 4 months ago by bk11 ★ 3.0k