Where communities thrive


  • Join over 1.5M+ people
  • Join over 100K+ communities
  • Free without limits
  • Create your own community
People
Repo info
Activity
    Paloma Contreras
    @PalomaUChile

    Hello, i'm trying to run the tutorial of selection tools and it isn't generating the file with iHS values. The problem is in this line of the pipeline:

    Rscript ../SelectionTools/selectionTools/corescripts/multicore_iHH.R -p CEU -i CEU_genetic_dist.haps -c 2 --window 5000000 --overlap 2000000 --maf 0.0 --big_gap 0 --small_gap 0 --small_gap_penalty 0 --haplo_hh --physical_map_haps CEU_genetic_dist.pos --cores 1 --working_dir . --offset 1 --ihs

    And the Error message is:
    Error in scan_hh(d, big_gap = opt$big_gap, small_gap = opt$small_gap, :
    unused arguments (big_gap = opt$big_gap, small_gap = opt$small_gap, small_gap_penalty = opt$small_gap_penalty, physical_positions = physical_positions)

    What can I do?
    murraycadzow
    @murraycadzow
    Hi,
    It looks like you need to install the version of rehh that comes with the pipeline because we modified it to the arguments for gap penalties
    Murray
    Paloma Contreras
    @PalomaUChile
    It worked, thanks for your answer!!
    Paloma Contreras
    @PalomaUChile
    Hi again,
    I have another question, but not as technical as the first one:
    Why the pipeline computes p-values, if the paper of iHS (Voith et al. 2006) says that it shouldn't be done?
    I couldn't find an explanation for this, not even in the rehh paper.
    Where can I find some information about this?
    Thanks for your help!
    murraycadzow
    @murraycadzow
    Hi,
    The pipeline produces p-values for iHS entirely because that is what the rehh package automatically does when standardising the iHH.
    As for if p-values should be calculated, this entirely depends on how you interpret them. My interpretation of why Voight et al 2006 doesn't like assigning p-values when they are used as the likelihood that a SNP has been selected which I too agree is not appropriate. However because of the relationship between a SNPs standardised iHS and its p-value it would be appropriate to use in a threshold manner ie P < 0.05 = the 5% most extreme iHS results
    Paloma Contreras
    @PalomaUChile

    Thanks for your answer! Sorry to bother you again, but now I'm having a problem running the pipeline.
    This time I'm using my own dataset and it worked well, except for three chromosomes (1, 9 and 16) that didn't complete the analysis correctly.
    The command that fails is:

    multcore_ihh command = Rscript /sfw/selectionTools/selectionTools/corescripts/multicore_iHH.R -p AYM -i AYM_genetic_dist.haps -c 9 --window 5000000
    --overlap 2000000 --maf 0.0 --big_gap 200000 --small_gap 20000 --small_gap_penalty 20000 --haplo_hh --physical_map_haps AYM_genetic_dist.pos --cores 40 --working_dir . --offset 1 --ihs

    once I run it, it returns this message:

    Error in [<-.data.frame(*tmp*, , 2, value = c(45440, 109810, 175050, :
    replacement has 14060 rows, data has 13391
    Calls: [<- -> [<-.data.frame

    I have checked the .haps and .pos file and they have the same amount of SNPs. Could it be a problem with some gaps in my dataset? I thought that the pipeline should just skip those regions (with gaps) and continue computing the iHS.

    I tried changing the window's size to 8Mb and it worked, but I'm not sure if that's the right way to solve it or if there is any better way to do it.
    I would really appreciate you advice in this matter. Thanks a lot!

    Regards,

    murraycadzow
    @murraycadzow

    Hi,
    I have had this problem in my own data and did change the way multicore_iHH.R rejoined the individual iHH files after the calculation step.
    The new method was updated into the selectionTools1.1 branch. It requires running the also updated haps_interpolate step before the updated multicore_iHH.R
    Changing window size is valid, ideally a larger window size is used where possible so that the ehh decay for the SNPs isn't prematurely stopped at the window boundary

    Murray

    Paloma Contreras
    @PalomaUChile
    I'm using the updated scripts now, and multicore_iHH.R doesn't work, not even with the chromosomes that were OK before.
    When I run this:
    Rscript /sfw/selectionTools/selectionTools/corescripts/multicore_iHH.R -p AYM -i AYM_genetic_dist.haps -c 16 --window 5000000 --overlap 2000000 --maf 0.0 --big_gap 200000 --small_gap 20000 --small_gap_penalty 20000 --haplo_hh --physical_map_haps AYM_genetic_dist.pos --cores 40 --working_dir . --offset 1 --ihs
    Paloma Contreras
    @PalomaUChile
    the error message is:
    Error in [.data.frame(map_positions, , 2) : undefined columns selected
    Calls: [ -> [.data.frame
    Execution halted
    (because the file AYM_genetic_dist.pos has 1 column even if I use the new script haps_interpolate to build it)
    Paloma Contreras
    @PalomaUChile
    I changed the Rscript so it can use the first column in the .pos file, and finally, I got the same problem as before:
    [1] "t_AYM.haps13.iHH does not exist! Continuing without"
    [1] "t_AYM.haps14.iHH does not exist! Continuing without"
    [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
    [26] 26 27 28 29 30 31 32
    Error in [<-.data.frame(*tmp*, , 3, value = 16) :
    new columns would leave holes after existing columns
    Calls: [<- -> [<-.data.frame
    Execution halted
    Paloma Contreras
    @PalomaUChile
    What can I do?
    James Boocock
    @theboocock

    This is to do with the markers missing names, something that was missed in the previous analyis.

    I can have a look today, will get back to you.

    James Boocock
    @theboocock
    Are you able to share a line from your haps file, this would enable me to reproduce the bug, making fixing it much easier.
    stevanspringer
    @stevanspringer
    Hi, I'm using selection tools on Docker on OSX. Thanks very much by the way, it's a fantastic tool. I've been mounting a folder called referencefiles with the -v flag. All neccesary subfolders and data pre-arranged on the local hard drive which allows me to get things in and out of the docker vm easily. It works very well.
    But, (there's always a but) the pipeline seems to be working great with some datasets, but it sometimes exits out analysis prematurely. Taking data from the 1000g for example. 6.89177428-89377427.
    stevanspringer
    @stevanspringer

    6.89177428-89377427 runs fine, but 19.20959639-21159638 fails, as does X.109745664-109945663. Test lactase dataset works well. selection_stderr.tmp is as follows:

    Traceback (most recent call last):
    File "/usr/local/bin/selection_pipeline", line 9, in <module>
    load_entry_point('selectionTools==1.0', 'console_scripts', 'selection_pipeline')()
    File "/usr/local/lib/python2.7/dist-packages/selectionTools-1.0-py2.7.egg/selection_pipeline/selection_pipeline.py", line 240, in main
    s.run_pipeline()
    File "/usr/local/lib/python2.7/dist-packages/selectionTools-1.0-py2.7.egg/selection_pipeline/standard_run.py", line 248, in run_pipeline
    fayandwus = self.variscan_fayandwus(haps2_haps)
    File "/usr/local/lib/python2.7/dist-packages/selectionTools-1.0-py2.7.egg/selection_pipeline/standard_run.py", line 536, in variscan_fayandwus
    start_position = int(start_pos)
    ValueError: invalid literal for int() with base 10: 'pos\n'

    Seems to be failing in variscan now, though earlier runs made it through to iHH and failed there. I'm initiating this with the CEU and YRI population files provided with the lactase test data, and using data downloaded directly from the 1000g browser.

    Sorry to bother you with what is certainly some mistake on my end. I'm just stumped since some coordinates make it through and others do not. Thank you! -Steve
    stevanspringer
    @stevanspringer
    I'm using defaults.cfg unmodified, and --phased-vcf. It
    It'll also run with --impute-split-size 5000, but fails or succeeds the same way.
    murraycadzow
    @murraycadzow
    Hi @stevanspringer , could you provide the command you used to get the 19.20959639-21159638 region please
    thanks
    Murray
    stevanspringer
    @stevanspringer

    Absolutely,

    http://browser.1000genomes.org/Homo_sapiens/Location/View?r=19:20959639-21159638
    click Get VCF data, confirm coordinates in text entry box, click next, right click to save and extract.

    snips from my actual data file:

    fileformat=VCFv4.1

    FILTER=<ID=PASS,Description="All filters passed">

    fileDate=20150218

    reference=ftp://ftp.1000genomes.ebi.ac.uk//vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz

    source=1000GenomesPhase3Pipeline

    snip ~250 lines. The following lines intentionally truncated to fit in the chat.

    CHROM POS ID REF ALT QUAL FILTER INFO FORMAT HG00096 HG00097

    19 20862396 esv3643896;esv3643897 A <CN0>,<CN2> 100 PASS AC=2,2;AF=0.000399361,0.000399361;AN=5008;CS=DUP_gs;END=20974015;NS=2504;SVTYPE=CNV;DP=13403;EAS_AF=0,0;AMR_AF=0,0.0029;AFR_AF=0,0;EUR_AF=0,0;SAS_AF=0.002,0;VT=SV GT 0|0 0|0
    19 20959646 rs73543372 C T 100 PASS AC=61;AF=0.0121805;AN=5008;NS=2504;DP=12049;EAS_AF=0;AMR_AF=0.0029;AFR_AF=0.0446;EUR_AF=0;SAS_AF=0;AA=.|||;VT=SNP GT 0|0 0|0
    19 20959677 rs541375182 T G 100 PASS AC=1;AF=0.000199681;AN=5008;NS=2504;DP=14387;EAS_AF=0;AMR_AF=0;AFR_AF=0.0008;EUR_AF=0;SAS_AF=0;AA=.|||;VT=SNP GT 0|0 0|0 0|0
    19 20959686 rs554447352 C G 100 PASS AC=1;AF=0.000199681;AN=5008;NS=2504;DP=14414;EAS_AF=0;AMR_AF=0;AFR_AF=0.0008;EUR_AF=0;SAS_AF=0;AA=.|||;VT=SNP GT 0|0 0|0 0|0

    Thanks very much!

    murraycadzow
    @murraycadzow
    thanks, looking into what the problem might be
    stevanspringer
    @stevanspringer
    Thank you! Great pipeline too, it's sure to be very useful to lots of folks.