Hi Seongmun,

I was using your script on several avian libraries and I found a couple of problems with the code in SEXCMD.R that you might want to check out:

Lines 118-127:

if (sextype == "XY"){
  ct <- data.frame(Length=do.call(rbind, marker_length_list[startsWith(x = names(marker_seq_list), prefix = "chrX")]), Count=do.call(rbind, marker_seq_list[startsWith(x = names(marker_seq_list), prefix = "chrX")]))
  ct <- cbind(ct, do.call(rbind, marker_seq_list[startsWith(x = names(marker_seq_list), prefix = "chrY")]))
  colnames(ct) <- c("Length", "Xcount", "Ycount")
}
if (sextype == "ZW"){
  ct <- data.frame(Length=do.call(rbind, marker_length_list[startsWith(x = names(marker_seq_list), prefix = "chrZ")]), Count=do.call(rbind, marker_seq_list[startsWith(x = names(marker_seq_list), prefix = "chrW")]))
  ct <- cbind(ct, do.call(rbind, marker_seq_list[startsWith(x = names(marker_seq_list), prefix = "chrW")]))
  colnames(ct) <- c("Length", "Zcount", "Wcount")
}

In the ZW system there is an inconsistency with respect to the XY system chunk of code. I think this line:
ct <- data.frame(Length=do.call(rbind, marker_length_list[startsWith(x = names(marker_seq_list), prefix = "chrZ")]), Count=do.call(rbind, marker_seq_list[startsWith(x = names(marker_seq_list), prefix = "chrW")]))
should be modified in:
ct <- data.frame(Length=do.call(rbind, marker_length_list[startsWith(x = names(marker_seq_list), prefix = "chrZ")]), Count=do.call(rbind, marker_seq_list[startsWith(x = names(marker_seq_list), prefix = "chrZ")]))
In this way it mimics the code for the XY system.

In lines 118, 123, 131, 140
if (sextype = "XY") and if (sextype = "ZW") should be replaced with if (sextype == "XY") and if (sextype == "ZW")

Thank you again for your tool :)

Valentina

Hello,
I was using your pipeline and came across an issue in the 1.extract_chr_from_REF.py script. It runs until line 28 but then in line 29 it gives me the following error:

Traceback (most recent call last):
File "PATH/SEXCMD/1.extract_chr_from_REF.py", line 29, in
for idx in xrange(0, len(chr_dict[CHR]), 50):
KeyError: 'chrX'

I was wondering to get some insight into why this might be going on. Thank you very much.

Hello,

I was wondering if it would be possible to have the sex marker files for human and chickens (the ones generated for the publication of this tool) within the Example files, it would be a good resource for people who want to use this tool on those organisms again.

Thanks for your help,
Valentina

I am trying to make my own marker file for humans (since the software does not provide the one used in the paper).

In step 1 I got the following error

Traceback (most recent call last):
File "../../../aDNA_programs/SEXCMD/1.extract_chr_from_REF.py", line 28, in <m
outfile.write('>' + CHR + '\n')
I am using python 2.7
I tried different reference genomes (from NCBI, 1000K, UCSN) and builds 37/38 but did not worked.
Is there a solution for this?

In order to move forward I downloaded the chromosomes X and Y from NCBI, and changed the names of the fasta to CHRX and CHRY (because I think the script in step 1 was also doing this).

Then, in step 2 I met an allocation limit problem with lastz:

FAILURE: in add_segment()
table size (4,869,542,152 for 101,448,794 segments) exceeds allocation limit of 4,294,967,279;
consider raising scoring threshold (--hspthresh or --exact) or breaking your target sequence into smaller pieces

At first I tried the lastz_32 (more RAM memory) instead of lastz. Gave me the same error and also it was too slow.

Then I added in the 2nd command the argument --hspthresh=4000 (3000 is the default)
It return the same error
I gradually (per 1000) increased the value and at 10000 it did not gave me an error but it stack for two days. It did the same with 20000.

May I have some help with these issues?

Thank you in advance!