isaacovercast
commented
3 years ago Ok, yeah this is unusual. I can see this in a pair-3rad dataset. Here is code to investigate this phenomenon in a vcf file:
invcf = "/path/to/amaranth_outfiles/amaranth.vcf"
vcf_df = pd.read_csv(invcf, skiprows=10, sep="\t")
# Walk through the df rows and populate a dictionary formatted like this:
# {locus_id:[CHROM, MIN_POS, MAX_POS]}
posdict2 = {}
cur_loc = ""
for idx in range(len(vcf_df)):
if not idx % 10000 and idx != 0:
print(idx, end=", ")
this_loc = vcf_df.iloc[idx]["ID"].split("_")[0]
if this_loc == cur_loc:
try:
next_loc = vcf_df.iloc[idx+1]["ID"].split("_")[0]
posdict2[cur_loc][2] = vcf_df.iloc[idx]["POS"]
except:
# Reached the end
posdict2[cur_loc][2] = vcf_df.iloc[idx]["POS"]
else:
cur_loc = this_loc
posdict2[cur_loc] = [vcf_df.iloc[idx]["#CHROM"], vcf_df.iloc[idx]["POS"], 0]
# Find consecutive loci that overlap
keys = list(posdict2.keys())
key_pairs = list(zip(keys, keys[1:]))
for k1, k2 in key_pairs:
# If consecutive loci are on different CHROM then continue
if posdict2[k1][0] != posdict2[k2][0]:
continue
# If loc1 max is greater than loc2 min, then they overlap
if posdict2[k1][2] >= posdict2[k2][1]:
print(k1, k2)This produces a list of consecutive loci that overlap in some fashion, looks like this:
loc155 loc156
loc171 loc172
loc485 loc486
loc847 loc848
loc864 loc865
loc952 loc953
loc959 loc960
loc967 loc968This is guaranteed to not be the smartest way to do this, but it works. I don't have any idea why this is happening yet, just wanted to document the way to figure out where it's happening.
I ran a reference aligned assembly in ipyrad v.0.9.68 with 150bp paired-end ddRAD data and got some partially overlapping loci returned as separate loci and SNPs on those loci are duplicated in the SNP data set.
Here are some examples:
Like the first two shown, there's usually (but not always) one long and one short locus and the short locus is completely or nearly overlapping with the long one. Out of 11476 returned loci, there are about 200 overlapping ones (i.e. 100 overlapping regions).
The same sample is never represented in both of the overlapping loci. Some samples will have data for one and some for the other. Here's an example from my VCF file, which shows that the overlapping portion is made up of only a few samples in one of the two loci. This is the pattern in all the cases I've looked at.
Is there a parameter I can adjust that will resolve these issues during the assembly instead of have to resolve them after the run, which is proving difficult?