about reference-panel in STITCH

185 views
Skip to first unread message

LH LV

unread,
Aug 10, 2018, 5:40:21 AM8/10/18
to STITCH imputation
Hi, Robbie. Sorry to bother you again. But I have some uncertainty about reference panel.
(1) how to specify the reference panel? Is it specified by --genfile=GENFILE?
(2) I found that value of --genfile is [default ""], it's empty. If I set the file by result from high depth sequencing and some addtional sites (Note this sites are all wild in high depth sequencing), will this operation have an effect on imputation?

Looking forward to your reply! Thank you.

Robbie Davies

unread,
Aug 13, 2018, 1:39:38 PM8/13/18
to LH LV, STITCH imputation
Hey,

1) Reference panel is specified by "reference_haplotype_file" and "reference_legend_file", and optionally "reference_sample_file" and "reference_populations". Have a look STITCH.R --help or ?STITCH in R, and the human examples in examples.R (which can be downloaded by running the file), and you should be able to figure out how to input the files and what their formats are (they should match one of the IMPUTE2 download packages)
Note that the reference panel is just used for the first iteration, so it helps speed up STITCH as you can usually run fewer iterations before an acceptable solution, but otherwise is not used. If you want to run using only the reference panel, that can be done (see the 3rd human example)

2) genfile only reports estimates of imputation accuracy and does otherwise not influence imputation accuracy (i.e. the genotypes are not used as part of the imputation process)

Best,
Robbie




--
You received this message because you are subscribed to the Google Groups "STITCH imputation" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stitch-imputat...@googlegroups.com.
To post to this group, send email to stitch-i...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/stitch-imputation/501b0f23-3311-403b-92c7-d7db3a0ea7ba%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
Message has been deleted

Robbie Davies

unread,
Aug 30, 2018, 10:35:18 PM8/30/18
to LH LV, STITCH imputation
GenomicConstructor is from SeqLib, not STITCH per-se. It's complaining it can't read a BAM file from 1:9837-2250350. Is chromosome 1 listed as "1" or "chr1" in the BAM? What STITCH uses has to match what the BAM uses

Best,
Robbie

On Thu, Aug 30, 2018 at 9:59 PM LH LV <laihu...@gmail.com> wrote:
Hi, Robbie. Sorry for replying so late. I have got the relational parameters("reference_haplotype_file" and "reference_legend_file").  And the helpful information about genfile you sent here is very useful for me and have changed my comprehension about genfile. Thank you.

Now I have got some error like before, but I have check this region and have numbers of SNP(about 7000).
[2018-08-31 09:18:36] Generate inputs
[2018-08-31 09:18:36] Error in get_sampleReadsRaw_from_SeqLib(useSoftClippedBases = useSoftClippedBases,  : 
  GenomicRegion constructor: Failed to set region for 1:9837-2250350
And here is my script:
Sys.setenv(PATH = paste0(getwd(),":",Sys.getenv("PATH")))
library(STITCH)
STITCH(
bamlist = "/results/list/s10w.bamlist",
outputdir = "/results/stitch_1000GP_2M",
method = "diploid",
originalRegionName = "1.1.2000000",
regenerateInput = TRUE,
regionStart = 1,
regionEnd = 2000000,
buffer = 250000,
chr = "1",
niterations = 10,
nCores = 24,
reference_haplotype_file = "/results/1000GP_Phase3/1000GP_Phase3_chr1.hap.gz",
reference_legend_file = "/results/1000GP_Phase3/1000GP_Phase3_chr1.legend.gz",
reference_sample_file = "/results/1000GP_Phase3/1000GP_Phase3.sample",
reference_populations = c("CHB","CHS","CDX"),
posfile = "/results/stitch_1000GP_2M/s10w.imputate.sites",
K = 10,
nGen = 8000,
refillIterations = NA)
Looking forward to your reply. Thanks.

在 2018年8月14日星期二 UTC+8上午1:39:38,Robert Davies写道:

LH LV

unread,
Aug 30, 2018, 11:56:40 PM8/30/18
to STITCH imputation
Hi, Robbie. The helpful information about you sent here is very useful for me and I got the normal operation. Thank you.

在 2018年8月31日星期五 UTC+8上午10:35:18,Robert Davies写道:
Reply all
Reply to author
Forward
0 new messages