from ipyrad.analysis.utils import Params
Hi @isaacovercast,
I'm trying to interpret a disparity between two parts of the *stats.txt file in the *_output directory produced by ipyrad v0.9.13. Here's a summary at the top:
total_filters applied_order retained_loci
total_prefiltered_loci 0 0 157483
filtered_by_rm_duplicates 8365 8365 149118
filtered_by_max_indels 10508 10508 138610
filtered_by_max_SNPs 10248 1621 136989
filtered_by_max_shared_het 262 19 136970
filtered_by_min_sample 123882 119546 17424
total_filtered_loci 153265 140059 17424And then in the locus coverage summary it looks like there's only 1125 loci:
locus_coverage sum_coverage
1 0 0
2 0 0
3 0 0
4 0 0
...
536 0 1125
537 0 1125Are there additional filters being applied to the data that are not summarized at the head of the *stats.txt file? I'd certainly take another 16000 loci if I could get them. Or is there something amiss in the summary at the top of the file, and I in fact have 1125 loci in the final dtaset? Thanks for any help!
------- ipyrad params file (v.0.9.13)-------------------------------------------
lig-denovo-c90 ## [0] [assembly_name]: Assembly name. Used to name output directories for assembly steps
/global/scratch/toconnor/radcap/combo ## [1] [project_dir]: Project dir (made in curdir if not present)
Merged: lig1-s12, lig2-s12, lig3-s12, lig4-s12, lig5-s12, lig6-s12 ## [2] [raw_fastq_path]: Location of raw non-demultiplexed fastq files
Merged: lig1-s12, lig2-s12, lig3-s12, lig4-s12, lig5-s12, lig6-s12 ## [3] [barcodes_path]: Location of barcodes file
Merged: lig1-s12, lig2-s12, lig3-s12, lig4-s12, lig5-s12, lig6-s12 ## [4] [sorted_fastq_path]: Location of demultiplexed/sorted fastq files
denovo ## [5] [assembly_method]: Assembly method (denovo, reference)
## [6] [reference_sequence]: Location of reference sequence file
pair3rad ## [7] [datatype]: Datatype (see docs): rad, gbs, ddrad, etc.
TAATT, ATGCA ## [8] [restriction_overhang]: Restriction overhang (cut1,) or (cut1, cut2)
5 ## [9] [max_low_qual_bases]: Max low quality base calls (Q<20) in a read
33 ## [10] [phred_Qscore_offset]: phred Q score offset (33 is default and very standard)
6 ## [11] [mindepth_statistical]: Min depth for statistical base calling
6 ## [12] [mindepth_majrule]: Min depth for majority-rule base calling
10000 ## [13] [maxdepth]: Max cluster depth within samples
0.9 ## [14] [clust_threshold]: Clustering threshold for de novo assembly
1 ## [15] [max_barcode_mismatch]: Max number of allowable mismatches in barcodes
2 ## [16] [filter_adapters]: Filter for adapters/primers (1 or 2=stricter)
80 ## [17] [filter_min_trim_len]: Min length of reads after adapter trim
2 ## [18] [max_alleles_consens]: Max alleles per site in consensus sequences
0.05 ## [19] [max_Ns_consens]: Max N's (uncalled bases) in consensus (R1, R2)
0.05 ## [20] [max_Hs_consens]: Max Hs (heterozygotes) in consensus (R1, R2)
100 ## [21] [min_samples_locus]: Min # samples per locus for output
0.2 ## [22] [max_SNPs_locus]: Max # SNPs per locus (R1, R2)
6 ## [23] [max_Indels_locus]: Max # of indels per locus (R1, R2)
0.5 ## [24] [max_shared_Hs_locus]: Max # heterozygous sites per locus
0, 0, 0, 0 ## [25] [trim_reads]: Trim raw read edges (R1>, <R1, R2>, <R2) (see docs)
0, 0, 0, 0 ## [26] [trim_loci]: Trim locus edges (see docs) (R1>, <R1, R2>, <R2)
l, p, s, n, k, a, g, G, u, v, t, m ## [27] [output_formats]: Output formats (see docs)
## [28] [pop_assign_file]: Path to population assignment file
## [29] [reference_as_filter]: Reads mapped to this reference are removed in step 3
Hi @isaacovercast -- some more followup.
I was scanning through the .loci file and noticed that although I set max_Indels_locus to 6, it looked like I had at most 1 or 2 (and usually no) indels in the loci that were written out.
So I increased max_Indels_locus to 20 and found that the header of the stats file now agreed with information reported farther down:
total_filters applied_order retained_loci
total_prefiltered_loci 0 0 157483
filtered_by_rm_duplicates 9071 9071 148412
filtered_by_max_indels 718 718 147694
filtered_by_max_SNPs 10997 10364 137330
filtered_by_max_shared_het 287 34 137296
filtered_by_min_sample 139008 134314 2982
total_filtered_loci 160081 154501 2982 locus_coverage sum_coverage
1 0 0
2 0 0
...
534 0 2982
535 0 2982
536 0 2982
537 0 2982Confusingly, when I changed max_Indels_locus back 6, I get the same behavior I mentioned before: many more loci reported in the header than are in the final assembly. If I bump max_Indels_locus to 100 I get more loci, but some of them look pretty fishy (as expected).
Not sure if that helps -- perhaps I am misinterpeting or misusing the max_indels_locus parameter.
Hi all, a couple questions:
1) I am aligning my RAD data to a reference genome that is divergent from my study group and using the ‘denovo+reference’ option. I noticed that none of the output files contain the sequences or snp calls for the reference genome. Is there an option to have ipyrad include that information in the output files? Or do you have any suggestions about the easiest way for me to get that information out of ipyrad?
2) With respect to the vcf file output format, what are the variant calls in the ‘REF’ column? I'm not sure they are from the reference genome because they always seem to be present in the file, whether the locus is aligned to a reference or not. Also, in the ‘CHROM’ column of the vcf, all the variants are assigned a locus "number”. Is there an easy way to figure out the association between these loci numbers and scaffold number they map to on the reference genome?
Thanks for your help!