text:
CNVKit does not output all the accessible regions in the targets bed file
This question was also asked on Biostars
I am using CNVkit on my data using hg38 as reference. The command that I am using is the following:
cnvkit.py batch sample.bam  -n control.bam -m wgs -f reference.fasta  --target-avg-size 1000 --output-dir results/

So, on the CNVkit documentation page is reported that

The batch --method wgs option uses the given reference genome’s sequencing-accessible regions (“access” BED) as the “targets” – these will be calculated on the fly if not provided.

Problem is, when it calculates accesibile regions from the reference, there are a lot of them listed in the log file that are not present in the final .bed output, and I cannot figure out why. Accessible regions are later divided in bins, and bins without enough coverage are discarded later in the process. These are the outputs I get on a short example:
Log file:
WGS protocol: recommend '--annotate' option (e.g. refFlat.txt) to help locate genes in output files.
NC_000001.10: Scanning for accessible regions
      Accessible region NC_000001.10:10000-177417 (size 167417)
      Accessible region NC_000001.10:227417-267719 (size 40302)
      Accessible region NC_000001.10:317719-471368 (size 153649)
      Accessible region NC_000001.10:521368-2634220 (size 2112852)
      Accessible region NC_000001.10:2684220-3845268 (size 1161048)
      Accessible region NC_000001.10:3995268-13052998 (size 9057730)
      Accessible region NC_000001.10:13102998-13219912 (size 116914)
      Accessible region NC_000001.10:13319912-13557162 (size 237250)
      Accessible region NC_000001.10:13607162-17125658 (size 3518496)
      Accessible region NC_000001.10:17175658-29878082 (size 12702424)
      Accessible region NC_000001.10:30028082-103863906 (size 73835824)
      Accessible region NC_000001.10:103913906-120697156 (size 16783250)
      Accessible region NC_000001.10:120747156-120936695 (size 189539)
      Accessible region NC_000001.10:121086695-121485434 (size 398739)
      Accessible region NC_000001.10:142535434-142731022 (size 195588)
      Accessible region NC_000001.10:142781022-142967761 (size 186739)
      Accessible region NC_000001.10:143117761-143292816 (size 175055)
      Accessible region NC_000001.10:143342816-143544525 (size 201709)
      Accessible region NC_000001.10:143644525-143771002 (size 126477)
      Accessible region NC_000001.10:143871002-144095783 (size 224781)
      Accessible region NC_000001.10:144145783-144224481 (size 78698)
      Accessible region NC_000001.10:144274481-144401744 (size 127263)
      Accessible region NC_000001.10:144451744-144622413 (size 170669)
      Accessible region NC_000001.10:144672413-144710724 (size 38311)
      Accessible region NC_000001.10:144810724-145833118 (size 1022394)
      Accessible region NC_000001.10:145883118-146164650 (size 281532)
      Accessible region NC_000001.10:146214650-146253299 (size 38649)
      Accessible region NC_000001.10:146303299-148026038 (size 1722739)
      Accessible region NC_000001.10:148176038-148361358 (size 185320)
      Accessible region NC_000001.10:148511358-148684147 (size 172789)
      Accessible region NC_000001.10:148734147-148954460 (size 220313)
      Accessible region NC_000001.10:149004460-149459645 (size 455185)
      Accessible region NC_000001.10:149509645-205922707 (size 56413062)
      Accessible region NC_000001.10:206072707-206332221 (size 259514)
      Accessible region NC_000001.10:206482221-223747846 (size 17265625)
      Accessible region NC_000001.10:223797846-235192211 (size 11394365)
      Accessible region NC_000001.10:235242211-248908210 (size 13665999)
      Accessible region NC_000001.10:249058210-249240621 (size 182411)
NT_113878.1: Scanning for accessible regions
      Accessible region NT_113878.1:0-106433 (size 106433)
NT_167207.1: Scanning for accessible regions
      Accessible region NT_167207.1:0-547496 (size 547496)
NC_000002.11: Scanning for accessible regions
      Accessible region NC_000002.11:10000-3529312 (size 3519312)
      Accessible region NC_000002.11:3579312-5018788 (size 1439476)
      Accessible region NC_000002.11:5118788-16279724 (size 11160936)
      Accessible region NC_000002.11:16329724-21153113 (size 4823389)
      Accessible region NC_000002.11:21178113-31705550 (size 10527437)
      Accessible region NC_000002.11:31705551-31725939 (size 20388)
      Accessible region NC_000002.11:31726790-31816827 (size 90037)
      Accessible region NC_000002.11:31816828-31816854 (size 26)
      Accessible region NC_000002.11:31816855-31816858 (size 3)
      Accessible region NC_000002.11:31816859-33092197 (size 1275338)
      Accessible region NC_000002.11:33093197-33141692 (size 48495)
      Accessible region NC_000002.11:33142692-87668206 (size 54525514)
      Accessible region NC_000002.11:87718206-89630436 (size 1912230)
      Accessible region NC_000002.11:89830436-90321525 (size 491089)
      Accessible region NC_000002.11:90371525-90545103 (size 173578)
      Accessible region NC_000002.11:91595103-92326171 (size 731068)
      Accessible region NC_000002.11:95326171-110109337 (size 14783166)
      Accessible region NC_000002.11:110251337-149690582 (size 39439245)
      Accessible region NC_000002.11:149790582-234003741 (size 84213159)
      Accessible region NC_000002.11:234053741-239801978 (size 5748237)
      Accessible region NC_000002.11:239831978-240784132 (size 952154)
      Accessible region NC_000002.11:240809132-243102476 (size 2293344)
      Accessible region NC_000002.11:243152476-243189373 (size 36897)
NC_000003.11: Scanning for accessible regions
      Accessible region NC_000003.11:60000-17137943 (size 17077943)
NT_113878.1: Joining over small gaps
NT_167207.1: Joining over small gaps
Wrote GCF_000001405.25_GRCh37.p13_genomic.bed with 2 regions

So, multiple accessible regions are present for sequences NC_000001.10 , NC_000002.11 and NC_000003.11, but only regions for NT_113878.1 and NT_167207.1 are reported:
NT_113878.1     0       106433
NT_167207.1     0       547496

Why is that so? I cannot really figure out.
I am sure the missed regions have data, because I am also using another CNVs caller on the same exact data, and it shows high coverage for the same regions CNVkit is discarding. Also the coverage for control data is high in those discarder regions, so this shouldn't be the matter, but I have no ideas at this point.
EDIT:
The question was answered on Biostars. Looking into the code, CNVkit is skipping chromosomes which don't have canonical names, like 'NCs' (so, the fact that it outputs the 'NTs' chromosomes is kind of incorrect, because they're also non-canonical and they should be skipped).

Did you check if the regions for which your bam file has no coverage/low mapping score actually have data? I don't know this, but it seems at least possible that it will simply not report regions that have no data and then filter out regions that have low coverage/bad score. And yes, I'm clutching at straws here.

I haven't used CNVKit, but does it rely on having a set of control samples? Assuming it does, like most non-WGS CNV callers, could it be that it is excluding regions for which the controls don't have good data?

Oh hey, thank you for taking the time to update with the answer! It would be really, really great if you could post that as an actual answer so this can be marked as answered. It will also give you a reputation boost, making more of the site features available to you.

The OP received an answer on their question on Biostars:

The question was answered on Biostars. Looking into the code, CNVkit is skipping chromosomes which don't have canonical names, like 'NCs' (so, the fact that it outputs the 'NTs' chromosomes is kind of incorrect, because they're also non-canonical and they should be skipped).


text:
Feeding the output of one task into another in WDL
Hi there I'm new to WDL and I need to complete a two-step workflow,
Basically, what I need to do is to extract and merge two reads dataset (forward and reverse) for a sample and pass the merged, uncompressed file to a different tool for genome inference. I'm able to complete the first step with no major issues; however, when it comes to get my output to the next task I'm facing some problems. Below a copy of the workflow I wrote
version 1.0

workflow step2 {
    input {
        String PANGENIE_CONTAINER = "overcraft90/eblerjana_pangenie:2.0.1"

        File FORWARD_FASTQ # compressed R1
        File REVERSE_FASTQ # compressed R2
        String NAME = "sample" # how to loop over samples' name in numerical order (maybe grub names' prefix)!?

        File PANGENOME_VCF # input vcf with variants to be genotyped
        File REF_GENOME # reference for variant calling
        File FASTQ_FILE # sample's FASTQ
        String VCF_PREFIX = "genotype" # string to attach to a sample's genotype
        String EXE_PATH = "/app/pangenie/build/src/PanGenie" # path to PanGenie executable in Docker

        Int CORES = 24 # number of cores to allocate for PanGenie execution
        Int DISK = 300 # storage memory for output files
        Int MEM = 100 # RAM memory allocated
    }

    call reads_extraction_and_merging {
        input:
        in_container_pangenie=PANGENIE_CONTAINER,
        in_forward_fastq=FORWARD_FASTQ,
        in_reverse_fastq=REVERSE_FASTQ,
        in_label=NAME, #later can be plural
        in_cores=CORES,
        in_disk=DISK,
        in_mem=MEM
    }

    call genome_inference {
        input:
        in_container_pangenie=PANGENIE_CONTAINER, # not sure whether Docker needs to be re-run
        in_pangenome_vcf=PANGENOME_VCF,
        in_reference_genome=REF_GENOME,
        in_executable=EXE_PATH,
        fastq_file=FASTQ_FILE,
        prefix_vcf=VCF_PREFIX,
        in_cores=CORES,
        in_disk=DISK,
        in_mem=MEM
    }

    output {
        File sample = reads_extraction_and_merging.fastq_file
        File genotype = genome_inference.vcf_file
    }
}

task reads_extraction_and_merging {
    input {
        String in_container_pangenie
        File in_forward_fastq
        File in_reverse_fastq
        String in_label
        Int in_cores
        Int in_disk
        Int in_mem
    }
    command <<<
        cat ~{in_forward_fastq} ~{in_reverse_fastq} | gzip -dc > ~{in_label}.fastq
        cp ~{in_label}.fastq /home/mat/cromwell-executions/step2/
    >>>
    output {
        File fastq_file = "~{in_label}.fastq"
    }
    runtime {
        docker: in_container_pangenie
        memory: in_mem + " GB"
        cpu: in_cores
        disks: "local-disk " + in_disk + " SSD"
    }
}

task genome_inference {
    input {
        String in_container_pangenie
        File in_reference_genome
        File in_pangenome_vcf
        String in_executable
        String prefix_vcf
        File fastq_file
        Int in_cores
        Int in_disk
        Int in_mem
    }
    command <<<
        echo "vcf: ~{in_pangenome_vcf}" > /app/pangenie/pipelines/run-from-callset/config.yaml
        echo "reference: ~{in_reference_genome}" >> /app/pangenie/pipelines/run-from-callset/config.yaml
        echo $'reads:\n sample: ~{fastq_file}' >> /app/pangenie/pipelines/run-from-callset/config.yaml
        echo "pangenie: ~{in_executable}" >> /app/pangenie/pipelines/run-from-callset/config.yaml
        echo "outdir: /app/pangenie" >> /app/pangenie/pipelines/run-from-callset/config.yaml
        cd /app/pangenie/pipelines/run-from-callset
        snakemake --cores ~{in_cores}
    >>>
    output {
        File vcf_file = "~{prefix_vcf}.vcf"
    }
    runtime {
        docker: in_container_pangenie
        memory: in_mem + " GB"
        cpu: in_cores
        disks: "local-disk " + in_disk + " SSD"
        preemptible: 1 # can be useful for tools which execute sequential steps in a pipeline generating intermediate outputs
    }
}

I try to summarize what's going on. After getting my two reads dataset, merging them and unzipping the output to ~{in_label}.fastq; I have to get this new file and feed it to the next task, so my idea was to copy it to a directory and form there get its path in the .json I use to run the WDL with Cromwell. First of all, I believe it might not be the right way to do it and if it is the case what should I do instead? Secondly, is there anything else in the code which might be changed? For instance, I realised I cannot use bgzip to speed up the unzipping step taking advantage of multiple threads; this will be useful when scaling to multiple samples.
Thanks in advance!

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
File I/O error using nglview.show_biopython(structure)
So I have been trying to get into visualizing proteins in python, so after some research I ended up on a tutorial that was teaching you how to visualize a protein from the COVID-19 virus, so I went and setup anaconda, got jupyter notebook working in vscode, and downloaded the necessary files from the PDB database, and made sure they were in the same directory as my notebook but when I run the the nglview.show_biopython(structure) function I get ValueError: I/O opertation on a closed file. I'm stymied this is my first time using jupyter notebook so maybe there is something I'm missing, I don't know.
This what the code looks like
from Bio.PDB import *
import nglview as nv

parser = PDBParser()
structure = parser.get_structure("6YYT", "6YYT.pdb")
view = nv.show_biopython(structure)

The error looks like this
Output exceeds the size limit. Open the full output data in a text editor
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
~\AppData\Local\Temp\ipykernel_8776\2743687014.py in <module>
----> 1 view = nv.show_biopython(structure)

c:\Users\jerem\anaconda3\lib\site-packages\nglview\show.py in show_biopython(entity, **kwargs)
    450     '''
    451     entity = BiopythonStructure(entity)
--> 452     return NGLWidget(entity, **kwargs)
    453
    454

c:\Users\jerem\anaconda3\lib\site-packages\nglview\widget.py in __init__(self, structure, representations, parameters, **kwargs)
    243         else:
    244             if structure is not None:
--> 245                 self.add_structure(structure, **kwargs)
    246
    247         if representations:

c:\Users\jerem\anaconda3\lib\site-packages\nglview\widget.py in add_structure(self, structure, **kwargs)
   1111         if not isinstance(structure, Structure):
   1112             raise ValueError(f'{structure} is not an instance of Structure')
-> 1113         self._load_data(structure, **kwargs)
   1114         self._ngl_component_ids.append(structure.id)
   1115         if self.n_components > 1:
...
--> 200         return io_str.getvalue()
    201
    202

ValueError: I/O operation on closed file

I only get this error when using nglview.show_biopython, when I run the get_structure() function it appears to read the file just fine. I can visualize other molucles just fine, or maybe that's because I was using the ASE library instead of a file. I don't know, that's why I'm here.
Update: Recently I found out that I can visualize the protein using nglview.show_file() instead of using nglview.show_biopython(). Even though I can visualize proteins now and technically my problem has been solved I would still like to know why the show_biopython() function isn't working properly.

I tried this in Jupyter Lab and VS code but couldn't reproduce the error. What OS are you on and what are your versions of Python, Biopython and nglview? I would like to help fix this as I added the show_biopython function many moons ago. In particular it seems that the StringIO object at https://github.com/nglviewer/nglview/blob/master/nglview/adaptor.py#L193-L200 gets closed early, preventing it being read.

@jgreener I'm on windows, the anaconda python kernel version is 3.9.13, my nglview version is 3.03 and my biopython version is 1.80.

answered, better tried to here : https://stackoverflow.com/questions/74737766/file-i-o-error-using-nglview-show-biopythonstructure/74782007#74782007

I also figured out another way to fix this problem. After going back to the tutorial I was talking about I saw that it was made back in 2021. After seeing this I wonder if we were using the same verions of each package, turns out we were not. I'm not sure what version of nglview they were using, but they were using biopython 1.79 which was the latest verion back in 2021 and I was using biopython 1.80. While using biopython 1.80 I was getting the error seen above. But now that I'm using biopython 1.79 I get this

So I guess there is something going on with biopython 1.80, so I'm going to stick with 1.79

Thanks, for posting back a solution. Could you kindly upvote @pippo1980's response please? I'm not an expert, by any means, but it looks to be a good response to me. Certainly useful.

yep, same happens with my error, it goes away if I use biopython 1.79

I had a similar problem with:
from Bio.PDB import *
import nglview as nv

parser = PDBParser(QUIET = True)
structure = parser.get_structure("2ms2", "2ms2.pdb")

save_pdb = PDBIO()
save_pdb.set_structure(structure)
save_pdb.save('pdb_out.pdb')

view = nv.show_biopython(structure)
view

error was like in question:
.................site-packages/nglview/adaptor.py:201, in BiopythonStructure.get_structure_string(self)
    199 io_str = StringIO()
    200 io_pdb.save(io_str)
--> 201 return io_str.getvalue()

ValueError: I/O operation on closed file

I modified site-packages/nglview/adaptor.py:201, in  BiopythonStructure.get_structure_string(self):
def get_structure_string(self):
        from Bio.PDB import PDBIO
        from io import StringIO
        io_pdb = PDBIO()
        io_pdb.set_structure(self._entity)
        io_str = StringIO()
        io_pdb.save(io_str)
        return io_str.getvalue()

with :
def get_structure_string(self):
        from Bio.PDB import PDBIO

        import mmap

        io_pdb = PDBIO()

        io_pdb.set_structure(self._entity)

        mo = mmap_str()

        io_pdb.save(mo)

        return mo.read()

and added this new class mmap_str() , in same file:
import mmap
import copy

class mmap_str():

    import mmap #added import at top

    instance = None

    def __init__(self):

        self.mm = mmap.mmap(-1, 2)

        self.a = ''

        b = '\n'

        self.mm.write(b.encode(encoding = 'utf-8'))

        self.mm.seek(0)

        #print('self.mm.read().decode() ',self.mm.read().decode(encoding = 'utf-8'))

        self.mm.seek(0)

    def __new__(cls, *args, **kwargs):
        if not isinstance(cls.instance, cls):
            cls.instance = object.__new__(cls)
        return cls.instance

    def write(self, string):

        self.a = str(copy.deepcopy(self.mm.read().decode(encoding = 'utf-8'))).lstrip('\n')

        self.mm.seek(0)

        #print('a -> ', self.a)

        len_a = len(self.a)

        self.mm = mmap.mmap(-1, len(self.a)+len(string))

        #print('a :', self.a)

        #print('len self.mm ', len(self.mm))

        #print('lenght string : ', len(string))

        #print(bytes((self.a+string).encode()))

        self.mm.write(bytes((self.a+string).encode()))

        self.mm.seek(0)

        #print('written once ')

        #self.mm.seek(0)

    def read(self):

        self.mm.seek(0)

        a = self.mm.read().decode().lstrip('\n')

        self.mm.seek(0)

        return a

    def __enter__(self):

        return self

    def __exit__(self, *args):

        pass

if I uncomment the print statements I'll get the :
IOPub data rate exceeded.
The notebook server will temporarily stop sending output
to the client in order to avoid crashing it.

error , but commenting them out I get:

while using thenglview.show_file(filename) I get:

tha's because, as could be seen looking at the pdb_out.pdb file
outputted by my code, Biopytho.PDB.PDBParser.get_structure(name , filename) doesnt retrieve the pdb header responsible for generate full CRYSTALLOGRAPHIC SYMMETRY/or biopython can't handle it (not sure about this, help if you know better),  but just the coordinates.
Still don't understand what is going on with the :
--> 201 return io_str.getvalue()

ValueError: I/O operation on closed file

it could be something related to jupiter ipykernal ? hope somebody could shed more light into this, dont know how the framework runs, but is definitely different from a normal python interpreter. As an example:

Same code in one of my Python virtualenv, will run forever, so it could be ipykernel dont like StringIO()s or do something strange to them ?
OK thanks to the hint in the answer below, I went inspecting PDBIO.py in github repo for version Biopython 1.80 and compared the save method of PDBIO : def save(self, file, select=_select, write_end=True, preserve_atom_numbering=False): with the one in Biopython 1.79,
see first bit:

and last bit:

so apparently the big difference is the with fhandle: block in version 1.80.
So I realized that changing adaptor.py with adding a subclass of StringIO that looks like:
from io import StringIO
class StringIO(StringIO):

    def __exit__(self, *args, **kwargs):

        print('exiting from subclassed StringIO !!!!!')

        pass

and modifying def get_structure_string(self): like this:
def get_structure_string(self):
        from Bio.PDB import PDBIO
        #from io import StringIO
        io_pdb = PDBIO()
        io_pdb.set_structure(self._entity)
        io_str = StringIO()
        io_pdb.save(io_str)
        return io_str.getvalue()

was enough to get my Biopython 1.80 work in jupiter with nglview.
That told I am not sure what are the pitfalls of not closing the StringIO object we use for the visualization, but apparently its what Biopython 1.79 was doing like my first answer using a mmap object was doing too (not closing the mmap_str)

What operating system are you on? I'm wondering if it could be a Windows issue, I can also test later on my Windows machine. There was a similar error reported a while ago for iPython and StringIO: https://github.com/ipython/ipython/issues/9168.

I am on a Debian VM (VMware on Debian Host), dont remember ipykernel version but tried downgrading it lower than 8 because of this  https://github.com/ipython/ipython/issues/9168 but getting same error nevertheless

This is an issue with Biopython 1.80, in particular introduced by a commit that automatically closes IO objects sent to PDBIO.save.
I made a PR to change the behaviour which fixes this issue and should hopefully be in the next Biopython release. Until then, avoiding using Biopython 1.80 should work.
Thanks @jeremiah-wade and @pippo1980 for digging into this.

it took You 2 lines to explain what it took me a couple of days to figure out ... . My question is what is / if is there one the pitfall of having unclosed IO objects around ? I am really new to coding and python so in really plain words if you can. thanks again

and on the same line, any chance to modifying the nglview adaptor.py  class BiopythonStructure(Structure) ??

It's explained a bit in the linked PR but you can't read from a closed StringIO object. So if the PDBIO.save function closes the object then an error will be thrown when the string is read back. In general closing an IO object is best done with a with block or .close(), rather than being done silently when the object is given as a function argument. In general the with block is better since it's quite common to forget to call .close(). Closing IO objects once they are no longer needed frees memory and lets other software use them safely.

They could just add a if isinstance(filename , StringIO) and close_file = 0 in PDBIO like 1.79 and change close_file = 1 for the #filehandle I hope bit

think solution could be like PDBParser handles it: from Bio.File import as_handle ,  https://github.com/biopython/biopython/blob/7ef56557bdca66eabd62334e61c4dc3153c997ff/Bio/File.py#L29 and in PDBIO.py  save method change  with fhandle:  to with as_handle(fhandle):

or remove  if isinstance else ... and  change with fhandle:    to   with as_handle(file, mode = 'w') as fhandle:

tried to understand git, I ended up with this, seems more coherent with the previous habits in the biopython project, but cant push it.

Anyone could pass it along ? [of course needs to be checked out, my tests are OK, but I am a novice].

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
I have a partially described gene from a plant, which tool may I use for predicting regolatory and UTRs sequences from a single DNA fragment?
I have a partially described gene from a plant with CDS and 3'UTR described but no 5'UTR; if I know the locus and upstream and downstream regions, is there a tool I can use for predicting regulatory DNA and UTRs for a given gene?
My annotation data source is the assembly of Phaseolus vulgaris here.

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

I think you need to tell us whether the genome has been sequenced. If not the nearest species with a sequenced genome.

when you say "partially described" what do you mean? what's the name? Do you have any sources of e.g. GFF3 annotation files, a la M__'s question?

Hi, sorry, this is my first post here.  So, the plant is Phaseolus vulgaris and is sequenced yet. When I say partially described I mean that I don't have (or can't find) regulatory sequences and in the gff3 file there is not a 5' (the mRNA starts with the first codon). https://ftp.ensemblgenomes.org/pub/plants/release-55/gff3/phaseolus_vulgaris/                    Thank you for the support, I'm new to bioinformatics.

updating with notes from comments

Assuming that you know the gene name, you would be able to look it up in a transcriptome assembly presumably, or perhaps NCBI has a full length transcript (search txid3885[Organism:exp] with your gene name).
For nearby regulatory DNA, you could consult an ATAC-seq dataset or other similar, though I doubt that causal links of regulatory DNA to your gene will have been directly assayed.


text:
How to extract date of acquisiton from mzXML and mzML files with Python?
I would like to extract date of acquisiton from mzXML and mzML files, using Python. Is it possible e.g. with Pyteomics? pymzml is only for mzML files. Are there other libraries that can do that? I assume of course that the date of acquisition is recored in mzXML/mzML files.

Please [edit] your question and i) add small example input files (they need to be valid mzXML and mzML) and ii) the output you expect from that input so we can test any answers we come up with.

This question is really targeted at people who know what a mzXML or mzML file is and who use e.g. pyteomics. Also, I don't know how to attach a data file like this to a stackoverflow question. I think, only images can be uploaded. Thanks

You can simply paste a few lines into your question and then use the {} button to format as code. And it isn't a question of showing the format but of having clear input and output requirements so solutions can be tested.

I don't think that is necessary.

It isn't about necessary, it is about helping others help you. If you include a test file, then we can test our solutions and make sure they actually work. You have the option of making it harder for people to help you, so in that sense it isn't "necessary". But if you want to ask a good, well received question that doesn't get downvoted, then it is necessary. Your call.

You should be able to extract anything from mzML and mzXML with Pyteomics. You will have to look at your files and figure out what you need to do, though, because you need to know what element contains the information that you want. I don't often work with mzXML and in the files that I have I do not see an acquisition date at all. With mzML it appears that the <run> element has this information:
<run id="small_raw" defaultInstrumentConfigurationRef="IC1" startTimeStamp="2005-07-20T14:44:22">

So I can get it like this:
In [1]: from pyteomics import mzml

In [2]: with mzml.MzML('tests/test.mzML') as f:
   ...:     print(next(f.iterfind('run', recursive=False))['startTimeStamp'])
   ...:
2005-07-20T14:44:22

P.S. Disclaimer: I am the maintainer of Pyteomics.

pymzml certainly works via Bald et al 2012:
conda install -c bioconda pymzml

Whilst this only works with mzML you simply convert mzXML to pymzml via ProteoWizard
msconvert your_mzXML_file --mzML

It is described here with your previous question with @haci
I believe you managed this via
conda create -n openms -c bioconda openms
FileConverter -in file.mzXML -out file.mzML

I suspect the underlying question is direct parsing in mzXML without the need to perform a conversion from mzXML to mzML or without using Bioconductor

Possibility The following may be of use massspec-toolbox at:
http://code.google.com/p/massspec-toolbox/source/browse/#svn/trunk/mzxml
Note, I've not tested nor read this but the repository looks genuine.


text:
How do I install shapeit4 on Ubuntu/Mint?
I'm having trouble install shapeit4 4.2.2 on my new Mint 21 installation (Ubuntu 22.04 LTS). According to the official docs, the requirements are:

HTSlib: A great C library for reading/writing high-throughput sequencing data.
BOOST: A free peer-reviewed portable C++ source libraries. SHAPEIT4 uses two specific BOOST libraries: iostreams and program_options.

I installed HTSlib with make and sudo make install, and libboost using its calls ./bootstrap and ./b2 but I still get errors when trying to run the make file for shapeit4:
/usr/bin/ld: cannot find /home/user/Tools/htslib-1.11/libhts.a: No such file or directory
/usr/bin/ld: cannot find /usr/lib/x86_64-linux-gnu/libboost_iostreams.a: No such file or directory
/usr/bin/ld: cannot find /usr/lib/x86_64-linux-gnu/libboost_program_options.a: No such file or directory

Manually editing makefiles is just... ugh. Save yourself the effort and use conda:
conda create -n shapeit
conda install -c bioconda shapeit4 -n shapeit

The documentation appears to be missing an edit one has to do. Open the makefile in the shapeit4 dir, and you will see something like (for version 4.2.2):
#COMPILER MODE C++11
CXX=g++ -std=c++11

#HTSLIB LIBRARY [SPECIFY YOUR OWN PATHS]
HTSLIB_INC=$(HOME)/Tools/htslib-1.11
HTSLIB_LIB=$(HOME)/Tools/htslib-1.11/libhts.a

#BOOST IOSTREAM & PROGRAM_OPTION LIBRARIES [SPECIFY YOUR OWN PATHS]
BOOST_INC=/usr/include
BOOST_LIB_IO=/usr/lib/x86_64-linux-gnu/libboost_iostreams.a
BOOST_LIB_PO=/usr/lib/x86_64-linux-gnu/libboost_program_options.a

So the paths appear to be hardcoded. Set these to the paths where you installed htslib (and compile it first).
For libboost, I found that its own installer doesn't place files in this expected location. So you can either change these paths to where you installed libboost, or you can use the apt manager version which will put the files in the right location by default. This is done with sudo apt-get install libboost-all-dev (found here).
After this, just run the make file again, and the output should looks like this:
user@computer:/dir/shapeit4-4.2.2$ make
make: Nothing to be done for 'all'.

Actually, precompiled SHAPEIT 4.2.2 is available in a Ubuntu repository.
sudo apt install shapeit4

Your answer could be improved with additional supporting information. Please [edit] to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center.


text:
Reading a Fast5-file with Python
I am trying to extract data from fast5-file with python 3.9.13 in Ubuntu. I have found a library "fast5_research"(This package comprises an API to HDF containers used by the research groups within Oxford Nanopore Technologies), installed on my PC and updated my software with "sudo apt update". According to tutorial I attempted to read fasta5-file with the next code:
from fast5_research import Fast5
filename='/home/master/Desktop/0.fast5'
with Fast5(filename) as fh:
    raw = fh.get_read(raw=True)
    summary = fh.summary()
print('Raw is {} samples long.'.format(len(raw)))
print('Summary {}.'.format(summary))

but as a result, I have an error:
KeyError                                  Traceback (most recent call last)
/tmp/ipykernel_410530/3650860803.py in <module>
      4 from fast5_research import Fast5
      5 filename='/home/master/Desktop/0.fast5'
----> 6 with Fast5(filename) as fh:
      7     raw = fh.get_read(raw=True)
      8     summary = fh.summary()

~/anaconda3/lib/python3.9/site-packages/fast5_research/fast5.py in __init__(self, fname, read)
     85
     86         # Attach channel_meta as attributes, slightly redundant
---> 87         for k, v in _clean_attrs(self[self.__channel_meta_path__].attrs).items():
     88             setattr(self, k, v)
     89         # Backward compat.

h5py/_objects.pyx in h5py._objects.with_phil.wrapper()

h5py/_objects.pyx in h5py._objects.with_phil.wrapper()

~/anaconda3/lib/python3.9/site-packages/h5py/_hl/group.py in __getitem__(self, name)
    175                 raise ValueError("Invalid HDF5 object reference")
    176         else:
--> 177             oid = h5o.open(self.id, self._e(name), lapl=self._lapl)
    178
    179         otype = h5i.get_type(oid)

h5py/_objects.pyx in h5py._objects.with_phil.wrapper()

h5py/_objects.pyx in h5py._objects.with_phil.wrapper()

h5py/h5o.pyx in h5py.h5o.open()

KeyError: 'Unable to open object (component not found)'

How to fix this error?

similar too https://stackoverflow.com/questions/28782487/read-hdf5-based-file-as-a-numpy-array-in-python

HDF5View should do the same as the post, just more "point and click"

The code looks fine fast5 is simply HDF5 format (rapid I/O).
The error,

'Unable to open object (component not found)'

You are not the only person to have this error and it could be a bug in the code with hdf5 in Python. Just to note I use hdf5 without a any problem.

Problem with the installation (?), this is the best guess to resolve it but I notice from your error log .. the errors are coming from Anaconda, I assume you made a conda installation, but the post implies apt-get?? Just check where the package is located.

conda install -c hcc fast5-research

Then perform a conda environment creation  and activation of the same package as here. I think you've already done that. anaconda3 is definitely your Python3 version ... and in that case it might be better to do the full conda installation because it stops things getting mixed up.
NOTE
There is a Docker container
docker pull quay.io/biocontainers/fast5-research

I would recommend using this (it will almost certainly work), if you know how to leverage Docker (which is sort of complicated).

File corruption ... examine the file via the HDFView GUI, here https://www.hdfgroup.org/downloads/hdfview/ . If you can't see within the file via the viewer - I think the file has a problem. Perhaps examine via hdf5 ... wait no that gets complicated.

It can't find the file or a path error (below).

Then kindly try:
from fast5_research import Fast5
import pathlib as Path
import sys

filename=Path('/home/master/Desktop/0.fast5')
if filename.exists():
    print ('Okay {} is fine.'.format(filename))
    filename='/home/master/Desktop/0.fast5'
else:
    pathway = Path(__file__).absolute().parent
    if not Path(pathway, '0.fast5').is_file():
        raise ValueError ('Sorry filename has not been found')
        sys.exit()
    filename = str(Path(pathway, '0.fast5'))
try:
    if Path(filename).exists():
        with Fast5(filename) as fh:
            raw = fh.get_read(raw=True)
            summary = fh.summary()
        print('Raw is {} samples long.'.format(len(raw)))
        print('Summary {}.'.format(summary))
except Exception as err:
    print ("The error codes are {} and {}".format(err.message, err.args))

Please execute the code in the same dir as '/home/master/Desktop/' or maybe just '/home/Desktop
In Fasta5 will almost certainly take a path from pathlib but I've erred in case not. HDF5 will certainly accept a pathlib pathway.
I've not tested it but just post back what happens.


text:
How do I determine read length from a user-inputted fastq file and then input that information into a tool's command line?
Can you guide me on how could take an input fastq file, read the length from read and then feed it into the findtail's command?
I have suggested concerns about documentation and program defaults to my prof and asked him to switch to cutadapt instead of using findtail that has no documentation. But he refuses so I have no option but to proceed with what he wants.
I want the user to input a RNA-seq data fastq file. Since the first line of each file contains the sequence length information, such as described:

@SRR3003130.1 HWI-ST741:596:H5FMVBCXX:1:1101:1244:1952 length=51

How do I pick the length=51, the 51 from this and feed in into the following command's --seqlength parameter:
    ./findtail --input_file Rrp6_S38_all_R1_001_egal.fastq --seqlength 51 \
--endgap 5 --taillength 10 --identity 100 --ptype A --stype T \
--output_format fasta --output_type pr > pr.fasta

The tool's command that I'm using only requires us to enter --seqlength once. So I do not think that determining the length of each and every read in the fastq file is of any use. Just simply the length from the first line of first read would suffice.
If something is lacking from my question, please let me know how can I improve.

And are you sure the length will be included in all files? What sequencers are your users using? Why not actually calculate the length? And what operating system are you using? The tools available depend on that.

I think that you cannot rely on that information being in a fastq header. Note also that each read has its own header! This is not just the first line! For example when I google "sample fastq file" the first examples do not have it. Spec doesn't seem to say anything about it.
For a one-off on a specific dataset this might be adequate:
grep -Eo "length=([0-9]+)" test.fq | cut -d"=" -f2
# prints "75" for each occurrence, so if you trust the first line you could pass to head:
grep -Eo "length=([0-9]+)" test.fq | cut -d"=" -f2 | head -1
# or you could explicitly sort and grab the longest read len:
grep -Eo "length=([0-9]+)" test.fq | cut -d"=" -f2 | sort | tail -1

For a more general application I'd suggest using some existing/explicit tools for measuring read length like seqkit or awk, rather than trusting some information that may or may not be in the header.
I can't find any docs for findtail, depending on what your application is I'd suggest using cutadapt instead, which does not require a single sequence length input.

To get the length of the first sequence from a single-line-per-sequence fastq, the second line can be piped through wc -c, and stored in a shell variable:
seqLength=$( head -n 2 input_file_name.fastq | tail -n 1 | wc -c)

This variable can then be used for the program arguments using ${variable}:
seqLength=$(head -n 2 input_file_name.fastq | tail -n 1 | wc -c)
./findtail --input_file input_file_name.fastq --seqlength ${seqLength} \
    --endgap 5 --taillength 10 --identity 100 --ptype A --stype T \
    --output_format fasta --output_type pr > pr.fasta

If the input file is gzipped [I have no idea if findtail supports gzipped files], then the length can be determined by piping through zcat:
seqLength=$( zcat input_file_name.fastq.gz | head -n 2 | tail -n 1 | wc -c)

As mentioned by many comments, note that some FASTQ files won't have the same length for each read, and it's unclear how findtail will behave for such a file.


text:
How to run how_are_we_stranded_here?
I'm trying to run the python library called how_are_we_stranded_here. I have installed with pip install how_are_we_stranded_here. I have also installed its dependencies i.e kallisto 0.44.0 and RSeQC. On how_are_we_stranded_here's pypi page, it gives the following command under the usage heading:
 check_strandedness --gtf Yeast.gtf --transcripts Yeast_cdna.fasta \
--reads_1 Sample_A_1.fq.gz --reads_2 Sample_A_2.fq.gz

So I tried running the program in the terminal using my own files, the command was:
    check_strandedness --gtf Glycine_max.Glycine_max_v2.1.55.gtf \
--transcripts Glycine_max.Glycine_max_v2.1.cdna.all.fa \
--reads_1 SRR2079644.fastq

But it keeps giving the error

./how_stranded.sh: line 1: check_strandedness: command not found

I will be very thankful for any help on how to run it.

You might need to include your pip install output too.  The package defines a set of scripts that should get you check_strandedness in your path, and it works for me... so I bet something didn't work right during install.

Would you consider installing this in a conda environment or other virtual environment, to control the background more effectively (with regard to M__'s answer and your comment below). It sounds like there are some issues in the background that might need resolving.

pip should just work, anyway... Two approaches, point 1 is a last resort approach (not recommended), I definitely recommend point 2.
First try,
python3 -m pip install --upgrade how_are_we_stranded_here

Try again does it work? Error logs?, then possibly ...

You need to find out where it is,

   pip show how_are_we_stranded_here
   pip show check_strandedness
   pip show how_stranded.sh

What you really want is the check_strandedness
cd ~
mkdir bin
./pathtoapp/check_stranded # make sure it works
ln -si /pathwaytoapp/how_stranded.sh /home/whatever/bin/
ln -si /pathwaytoapp/check_strandedness ~/bin/check_strandedness

Then write
export PATH="/home/whatever/bin:$PATH"

or modify .profile or .bashrc, then
source ~/.bash_profile

Try again ...

It is better to perform a conda installation

Install either,

Miniconda or;
Anaconda

    conda install -c bioconda how_are_we_stranded_here
    conda create --name stranded
    conda activate stranded
    check_strandedness -h

Thank you so much for taking out the time to give such a thorough reply. It definitely solved the check-strandedness issue but opened a plethora of other issues, so I just gave up.

No problem @pubsurfted, although this one is not 'solved' upvotes always welcome.

upvoted. Thank you!


text:
How to remove double '/' from file path, bash script
I'm using the following script to detect strandedness of my paired end RNA-seq data.
#!/usr/bin/env bash

cd ~/Desktop/fastq

for infile in *_1.fastq
do
   base=$(basename ${infile} _1.fastq)

    docker run quay.io/biocontainers/guessmylt:0.2.5--py_0 GUESSmyLT \
    --reads ${infile} ${base}_2.fastq

done

But when I run it, it keeps giving the error:
Building DAG of jobs... File path //SRR3475890_1.fastq contains double
'/'. This is likely unintended. It can also lead to inconsistent
results of the file-matching approach used by Snakemake. File path
//SRR3475890_2.fastq contains double '/'. This is likely unintended.
It can also lead to inconsistent results of the file-matching approach
used by Snakemake. MissingInputException in line 46 of
/usr/local/lib/python3.6/site-packages/GUESSmyLT/subsample: Missing
input files for rule subsample_fastq: //SRR3475890_2.fastq
//SRR3475890_1.fastq

Why is it appending double '/' to the fastq file's path? Nowhere in my script did I mention any path.
I apologize for the very vague title, English is not my first language, so I did not know how else to phrase it.

The issue here is that the current working directory inside the container is /. You can see this with:
docker run --rm quay.io/biocontainers/guessmylt:0.2.5--py_0 python \
    -c 'import os; print(os.getcwd())'

When GUESSmyLT tries to get the absolute path of it's arguments (--reads, --reference, --output, etc), it actually appends a trailing slash to the current working directory1, hence the double slash in the returned file names.
To pass the files into the container and ensure that the results are available when it exits, you'll need a volume mount anyway. For example:
docker run --rm -v "${PWD}:/working" quay.io/biocontainers/guessmylt:0.2.5--py_0 GUESSmyLT \
    --reads "/working/${base}_1.fastq" "/working/${base}_2.fastq" \
    --reference "/working/reference.fa" \
    --annotation "/working/annotation.gtf" \
    --output "/working"

Ah of course, thanks. The citation to the GitHub Python code is cool.

Thank you @Steve for taking the time to give such a thorough reply. It works now!


text:
Grep peculiar pattern inside a file
I would like to grep everything inside the second square bracket of a file:
>AB1311658.1 unnamed protein product [C.sp. balc.]

>evm.model.pilon.8 [actin][D.r.]

>NP_0232343.1 [fucose mutarotase] [E.l.]

Obtaining:
C.sp. balc.

D.r.

E.l.

How can I solve it?
Marco

Is that everything you have in the file? No other lines? And what operating system are you using?

You can easily install GNU version of the commands into your Mac. Install the GNU Command Line Tools. You can do this using Homebrew, then
brew install coreutils. Default will be MacOS version (unless you intentionally set otherwise), but when you want to use GNU version, just add a "g" before the command. For example, "sort" vs "gsort"; "grep" vs "ggrep". It will make your life easier as a Mac user.

One way using sed:
sed -n '/^>/s/.*\[\([^]]*\)\]$/\1/p' input.fa

Results:
C.sp. balc.
D.r.
E.l.

Explanation:

Turn off automatic printing of the pattern space using -n
Address lines starting with a header using /^>/
Perform a substitution using a capture group to select everything between the last set of brackets on these lines, with s/.*\[\([^]]*\)\]$/\1/
Print the pattern space with p

You could probably simplify your sed to sed -n '/^>/{s/.*\[//; s/]//p}' file.

The /^>/{ doesn't appear to be necessary.

With GNU grep (the default on Linux):
$ grep -oP '.*\[\K[^]]+' file
C.sp. balc.
D.r.
E.l.

The -o means "print only the matched portion of the line" and -P enables Perl Compatible Regular Expressions which give us \K for "ignore everything matched up to this point". So we then look for the longest possible string ending in a [ and discard it (.*\[\K), and then look for the longest possible stretch of non-] characters after that.
For very large files, you can speed things up a little by focusing only on the header lines:
grep -oP '^>.*\[\K[^]]+' file

For OSX, I don't have a -P option (Perl pattern matching) so the best I can do is finish the one-liner by piping into perl
grep -o '\[.*\]\[[^]]*' input.fa | perl -p -e 's/[][]/ /g'

The -o makes grep match the substring, rather than the entire line.
With GNU grep (the default on Linux) and the -P option, it's a single line of grep where \K or (?=pattern), (?<=pattern) are used for non-capture characters and these are placed around the [
grep -oP '(?<=\[).*(?=\]).*(?<=\[)[^]]*' input.fa

The ?<= is a "positive lookbehind" and the ?= is a "positive lookahead".
It's a lot easier to use a perl one-liner.

You could also pipe to sed, and skip the Perl dependency: grep -o '\[[^]]*\]$' input.fa | sed 's/[][]//g'

@Steve what benefit would that give? You just added sed as a dependency instead. Both sed and perl are installed by default on the vast majority of systems, so unless you are using AIX or a very limited embedded system, perl is just as valid as sed.

@terdon I agree, but nearly all systems with grep will have sed.

@Steve yes, exactly like perl. Unless people are trying to do bioinformatics on toasters or something, I really don't see any point in using sed instead of perl. They are both very powerful, standard, ancient commands that have a long history in the *nix world in general and in bioinformatics in particular. The only full system (as opposed to embedded system) I know of that doesn't install perl by default is AIX and that isn't likely to be used for bioinformatics.

Perl can easily be downloaded @Steve ;), "Dockerised", or "conda-ised". Not many Linux OS have removed Perl, FreeBSD is one. Most of us work with Ubuntu, Debian or Darwin. I think even Alpine Linux has Perl.

@terdon Well, Perl is a sed-killer. For simple substitutions like the above, PCRE are not needed. You can use Perl if you want to, but you could also use sed and not need to escape the brackets.

@M__ Alpine does not come with Perl out of the box

@Steve well yes, but Alpine is mostly for embedded systems, right? Not sure what you mean about escaping, neither tool requires you to escape the brackets. And sure, of course you can use sed. I just don't see any benefit in using one external tool over another here. Nothing wrong with using sed, but nothing wrong with using perl either. Granted, if I bring out perl, I would probably do the whole thing in perl and remove the grep, but still.

@Steve yes technically Alpine Linux holds sed and grep but via Busybox and its grep is basic without -P. I installed perl on Alpine it took a few seconds  maximum. However, the reason is the container is 7.5 MG and Alpine Perl is 55MB. The key to the question is GNU grep.

@terdon You're right, neither tool requires the brackets to be escaped but 's/[[]]/ /g' should be 's/[][]/ /g I think.

@Steve eeeek! You're quite right, thanks!

@Steve here's an article specially written for you ;) https://stackoverflow.blog/2022/07/06/why-perl-is-still-relevant-in-2022/

Thanks @M__ - I hadn't even heard of Perl 7 until now. But then I haven't really touched Perl for nearly a decade now. I just hope you're not confusing me for someone who dislikes Perl - that's certainly not the case. It's just that it's not needed here in my opinion. The Perl dependency also makes pipelining the solution tricky using the standard (i.e. non-mulled) quay.io biocontainers which provide the busybox suite. But that's not likely even a concern here..

@Steve "not confusing me for someone who dislikes Perl " ... just checking ;) Fun as this is the last 3 comments at least need to be deleted, probably more

@terdon Alpine is a "freak" Linux OS, just for Docker containers. I should have checked my information before introducing it here :)

Some of the solutions suggested above utilize GNU grep. If you are a Mac user, you can easily install GNU version of the commands into your Mac. You can do this by installing the GNU Command Line Tools. The easiest way is to use Homebrew, then brew install coreutils.
The default commands will still be MacOS version (unless you intentionally set otherwise), but when you want to use GNU version, you could just add a g before the command. For example, sort vs gsort, grep vs ggrep. You can also check manual of the command the same way, e.g. man gsort. Trust me this will make your life easier as a Mac user.
Also, in recent years default MacOS Terminal shell is zsh. If you desire you can switch it back to bash.

Thanks @Supertech perfect


text:
HTTP status code: 400 "Value must be an integer" error uploading JSON schema results to benchling
I am writing a benchling uploader script in python to upload data from a csv file in a JSON formatted schema to the benchling API platform. My function to create the JSON schema in the script looks like the following:
def create_antismash_schema(dataframe):
    """Create a json object for antismash output"""
    return {
        "assayResults": [
            {
                "fields": {
                    "entity": {"value": j['SBI ID']},
                    "sequence_id": {"value": j['sequence_id']},
                    "category": {"value": j['Category']},
                    "product": {"value": j['Product']},
                    "start": {"value": j['start']},
                    "end": {"value": j['end']},
                    "nucleotide_length": {"value": j['nucl_length']}
                },
                "schemaId": "assaysch_Z6B7TlWG",
                "projectId": 'src_tIUAcxDH'
            } for i, j in dataframe.iterrows()
        ]
    }

Now the data I have in benchling for the fields matches the keys in the schema and the values match the headers in the csv file I'm pulling the data from.
benchling results schema looks like:

csv file:

The error I am getting when I try to upload is the following:

Uploading Antismash data to Benchling...
  0%|                                                                                | 0/159 [00:00<?, ?it/s]{
    "assayResults": null,
    "errors": [
        {
            "fields": {
                "nucleotide_length": "Value must be an integer"
            },
            "global": []
        },
        {
            "fields": {
                "nucleotide_length": "Value must be an integer"
            },
            "global": []
        },
        {
            "fields": {
                "nucleotide_length": "Value must be an integer"
            },
            "global": []
        }
    ]
}
  0%|                                                                                | 0/159 [00:00<?, ?it/s]
Failed to upload summary.
HTTP status code: 400

The field is set to be an integer in benchling and all the values I'm pulling from in the csv file for the nucl_length column are integer values so not sure why this error is beig thrown. Any guidance here appreciated, thanks.

There is a major issue about the json format required, however just converting to json is easy because iterrows has to be pandas and pandas has its own output json format ... For example,
dataframe.to_json('benchling.json'', orient='records', lines=True)

If that fails then,
jsonOut = dataframe.to_json(orient='records')[1:-1].replace('},{', '} {')
with open('benchling.json', 'w') as fout:
    fout.write(jsonOut)

However, I don't think it would work and at a minimum I think the headers of the csv will need to change.

The error

"nucleotide_length": "Value must be an integer"

Looking at your code ....
"nucleotide_length": {"value": j['nucl_length']}

Sure this is right?
Why not try ...
"nucleotide_length": int(j['nucl_length'])

Thats what the error seems to be saying.


text:
How to include only the trimmed fastq files into a nextflow channel and exclude any secondary log files created?
I have the following simple nextflow script, that takes fastq paired reads and runs them through the atria adapter trimming.
// Declare syntax version
nextflow.enable.dsl=2

/*
 * pipeline input parameters
 */

params.reads = "$projectDir/data/*_{1,2}.fastq"
params.transcriptome_file = "$projectDir/data/Glycine_max.cdna.all.fasta"
params.gtf_file = "$projectDir/data/Glycine_max.gtf"
params.outdir = "results"

/*
 * define the `TRIM` process that creates trimmed.reads.fastq files
 * given the reads.fastq file.
 */

process TRIM {

    tag "$pair_id"

    input:
    tuple val(pair_id), path(reads)

    output:
    path(pair_id)

    script:
    """
    atria -r ${reads[0]} -R ${reads[1]} -t 1 -l 12 -q 30 -o $pair_id
    """
}

workflow {

    read_pairs_ch = Channel.fromFilePairs( params.reads, checkIfExists:true )

    trim_ch=TRIM(read_pairs_ch)

}

Next what I want to achieve is input the atria.trimmed.fastq files generated by the trim process into a channel so that I could use it as an input for other process. But atria trimming tool also creates two additional log files along with atria.trimmed.fastq. How do i exclude those additional files and only input the trimmed files into a channel for later use.
Best wishes.

Cross posted here: https://stackoverflow.com/questions/74780770/how-to-include-only-the-trimmed-fastq-files-into-a-nextflow-channel-and-exclude

With output: path(pair_id), you're telling Nextflow to expect a file or directory with the value of pair_id in the process working directory when the command exits (successfully). Atria outputs all files to a directory specified using -o PATH1. This directory is created if it doesn't already exist. Note that the output directory defaults to the current working directory if -o is not specified. If you'd like to write the results to a sub-directory, you could try:
process TRIM {

    tag { pair_id }

    publishDir "${params.outdir}/TRIM", mode: 'copy'

    cpus 1

    input:
    tuple val(pair_id), path(reads)

    output:
    tuple val(pair_id), path(trimmed_reads), emit: trimmed_fastqs
    path "${pair_id}/*.atria.log*", emit: logs

    script:
    def (r1, r2) = reads

    trimmed_reads = reads.collect {
        "${pair_id}/${it.baseName}.atria.fastq"
    }

    """
    atria \\
        -r "${r1}" \\
        -R "${r2}" \\
        -o "${pair_id}" \\
        -t ${task.cpus} \\
        -l 12 \\
        -q 30
    """
}

In the above, we used the emit option in the process output definition to assign a name identifier to each declaration. We can use this feature to reference the desired channel in the workflow block. For example:
workflow {

    read_pairs_ch = Channel.fromFilePairs( params.reads, checkIfExists:true )

    trim_out = TRIM(read_pairs_ch)

    trim_out.trimmed_fastqs.view()
}

Results:
$ nextflow run main.nf
N E X T F L O W  ~  version 22.10.0
Launching `main.nf` [distraught_dalembert] DSL2 - revision: aea2d37092
executor >  local (3)
[5f/78908b] process > TRIM (A) [100%] 3 of 3 ✔
[B, [/path/to/work/6a/440a50e9a0725c18f9eede80549cb8/B/B_1.atria.fastq, /path/to/work/6a/440a50e9a0725c18f9eede80549cb8/B/B_2.atria.fastq]]
[A, [/path/to/work/5f/78908b96adee8b8d663d9b4eb4faba/A/A_1.atria.fastq, /path/to/work/5f/78908b96adee8b8d663d9b4eb4faba/A/A_2.atria.fastq]]
[C, [/path/to/work/b0/a3bb3024a4c05fffec555cce5097a7/C/C_1.atria.fastq, /path/to/work/b0/a3bb3024a4c05fffec555cce5097a7/C/C_2.atria.fastq]]

To input the trimmed reads into a downstream process, just ensure that the input tuple matches the input set cardinality declared by the process. If the input declaration is not the same, you may need to apply one or more operators to your channel to make it fit. For your alignment process, you might not need to make any changes. Also, if your alignment process needs to declare one or more reference files, just make sure these are value channels. Most of the time, what you want is one queue channel and one or more value channels when you require multiple input channels, otherwise you might observe unusual behavior (e.g. see this question). The section on multiple input channels in the docs is well worth taking the time to read in my opinion.
A dummy alignment process might look like:
process ALIGN {

    tag { sample }

    debug true

    input:
    tuple val(sample), path(reads)
    path fasta
    path gtf

    script:
    def (r1, r2) = reads

    """
    ls -g "${r1}"
    ls -g "${r2}"
    ls -g "${fasta}"
    ls -g "${gtf}"
    """
}

And the updated workflow:
workflow {

    read_pairs_ch = Channel.fromFilePairs( params.reads, checkIfExists:true )

    transcriptome = file( params.transcriptome_file )
    gtf = file( params.gtf_file )

    TRIM( read_pairs_ch )

    ALIGN( TRIM.out.trimmed_fastqs, transcriptome, gtf )
}

Results:
$ nextflow run main.nf
N E X T F L O W  ~  version 22.10.0
Launching `main.nf` [extravagant_lattes] DSL2 - revision: 29b5315739
executor >  local (6)
[51/e4e261] process > TRIM (A)  [100%] 3 of 3 ✔
[cc/72c53d] process > ALIGN (A) [100%] 3 of 3 ✔
lrwxrwxrwx 1 users 57 Dec 14 00:21 C_1.atria.fastq -> ../../c1/9a85e19c601f2ab2c4af9c6fd203f0/C/C_1.atria.fastq
lrwxrwxrwx 1 users 57 Dec 14 00:21 C_2.atria.fastq -> ../../c1/9a85e19c601f2ab2c4af9c6fd203f0/C/C_2.atria.fastq
lrwxrwxrwx 1 users 40 Dec 14 00:21 Glycine_max.cdna.all.fasta -> ../../../data/Glycine_max.cdna.all.fasta
lrwxrwxrwx 1 users 29 Dec 14 00:21 Glycine_max.gtf -> ../../../data/Glycine_max.gtf

lrwxrwxrwx 1 users 57 Dec 14 00:21 B_1.atria.fastq -> ../../3f/16747fbd2c9bc4c0dc945b24873cca/B/B_1.atria.fastq
lrwxrwxrwx 1 users 57 Dec 14 00:21 B_2.atria.fastq -> ../../3f/16747fbd2c9bc4c0dc945b24873cca/B/B_2.atria.fastq
lrwxrwxrwx 1 users 40 Dec 14 00:21 Glycine_max.cdna.all.fasta -> ../../../data/Glycine_max.cdna.all.fasta
lrwxrwxrwx 1 users 29 Dec 14 00:21 Glycine_max.gtf -> ../../../data/Glycine_max.gtf

lrwxrwxrwx 1 users 57 Dec 14 00:21 A_1.atria.fastq -> ../../51/e4e261d019c65377d63826963b19b1/A/A_1.atria.fastq
lrwxrwxrwx 1 users 57 Dec 14 00:21 A_2.atria.fastq -> ../../51/e4e261d019c65377d63826963b19b1/A/A_2.atria.fastq
lrwxrwxrwx 1 users 40 Dec 14 00:21 Glycine_max.cdna.all.fasta -> ../../../data/Glycine_max.cdna.all.fasta
lrwxrwxrwx 1 users 29 Dec 14 00:21 Glycine_max.gtf -> ../../../data/Glycine_max.gtf

I appreciate your reply so much and I am very thankful for it. I wanted to confirm one more thing as I just started using nextflow. How do I input these trimmed_fastqs into a new process, say, for alignment? Do I use trim_out.trimmed_fastqs channel as input for the alignment process?

@pubsurfted Correct. I've tried to give an example above of what you might need. Please let me know if I've missed anything, thanks!


text:
Allele frequency gnomad data
I noticed that some gnomad SNPs have allele frequency =0 and some are just not found : No genes were found in this region. what is the difference between them : for example:Y-10010520-10010520.
Is all frq=0 means that it is very rare and therefore zero? then what about AF=NA? Also, does the gnomad genome data contain information about the exomes also, or is the information in the exome data is  not included in the genome data in gnomad?

Cross-posted on biostars: https://www.biostars.org/p/9548344/

I don't do eukaryotes, well humans anyway, but in population genetics an allele frequency of zero simply means it is present in no individuals against the reference allele. It doesn't mean the allele is heterozygous nor homozygous its simply absent. Whether that is informative or not is an interesting question. In addition, it does not necessarily mean the allele is rare, it is simply absent against the population in your analysis. The sickle cell anaemia allele is not rare, in certain parts of Africa it has high frequency, but it would be absent if all the individuals in the  population of the analysis that did not have African ancestry. In this specific example, it would likely reflect inaccurate sampling of the population in the analysis.
If allele frequency was 1 it would mean complete homozygosity and identity, which is very informative.
I'd simply refer to the wiki which is pretty good for this. Allele frequency is the entire basis of population genetics.


text:
How to keep the basename of an input file and change its extension in Nextflow?
The following Nextflow scripts takes the $gtf as an input and outputs a .bed file. What I want is to keep just the basename of the gtf file and add the .bed extension to it. I tried declaring a variable x and using basename from bash but it gives error.
// Declare syntax version
nextflow.enable.dsl=2

/*
 * pipeline input parameters
 */

params.gtf_file = "$projectDir/data/Glycine_max.gtf"

    /*
 * define the `GTF2BED` process that creates a BED file
 * given the params.gtf_file as input.
 */

process GTF2BED {

    input:
    path gtf

    output:
    "${gtf}/.bed"

    script:
    """
    gtf2bed --gtf $gtf --bed "${gtf}/.bed"
    """
}

Best wishes.

You can use the getBaseName() method on file objects to get the basename:

The invocation of any method name starting with the get prefix can be shortcut by omitting the get prefix and () parentheses. Therefore, writing myFile.getName() is exactly the same as myFile.name and myFile.getBaseName() is the same as myFile.baseName and so on.

For example:
process GTF2BED {

    input:
    path gtf

    output:
    path "${gtf.baseName}.bed"

    script:
    """
    gtf2bed --gtf "${gtf}" --bed "${gtf.baseName}.bed"
    """
}

If you need to remove a double extension (.e.g. .tar.gz) and don't want to remove all extensions (e.g. using the getSimpleName() method), you can use the long form to supply the number of extensions to remove. For example:
tarball = file( 'my.file.tar.gz' )

println( tarball.getBaseName(2) )

Prints:
my.file


text:
SNPs location out side the gene
I downloaded exome data from gnomad to get the allele frequency of some SNPs that I'm analyzing. The SNPs are from external resource. Some SNPs (in my data) are not present in the exome data (in gnomad):

"No genes were found in this region"

As part of my report question for homework I need to answer this:

"Are the SNPs in your data  that are not in gnomad are outside genes?",

What does it mean SNP outside genes? Does it mean that they are not in the exome data but in the genome data on gnomad?
I found this in an article:

"Single-nucleotide polymorphisms (SNPs) are single base-pair substitutions that occur within and outside genes"

Also found this:

"SNP markers might be selected outside of exomes or the direct cause of a genetic mutation" does this means that it is found in the genome but not exome data in gnomad as it is outside the gene?

Basically, I can't seem to understand what is a SNP outside a gene and can I see this information in gnomad.

@Community edited hope it is better

Cross-posted on biostars: https://www.biostars.org/p/9548287/

The answer is you can't see a SNP outside a gene for gnomad because you stated it is exome data. Thus, a minimal definition of the gene is it comprises: introns, exons and a promoter region.
Exome = exon so you downloaded a subsection the gene across the genome.
To reiterate, you don't have any external allele frequency data outside the gene. However, if your in-house data set is whole genome data, i.e. for making your SNP calls then most of these SNPs will be outside the exon (thats a result of comparative selection pressure). So the answer is yes. It isn't clear what the in-house data set comprises however and outside an exon does not necessarily mean its outside the gene. A lot of mutations will occur within introns (which is outside the exon but still within a gene), but in addition more mutation will occur in intergenic regions, i.e. completely outside the exon + intron + promoter. In this sense the question seems to possibly lack clarity.
To be honest, questions on bioinformatic practicals are often after-thoughts to the core activity, so I wouldn't worry too much about the semantics.

♦ thank you for the detailed answer. so if I understand correctly if instead of looking the exome data I would look at the whole genome data in gnomad that would be looking outside the gene?

Yes that would resolve the situation.


text:
Kallisto error: index input file could not be opened!
I am utilizing Kallisto in Anaconda/miniconda for RNA sequencing. I have successfully made an index; using said index to analyze RNA sequencing data has yielded an error of:

"Error: index input file could not be opened!"

The implication is that the code used does in fact work but raises the question as to why the index file could not be opened. I have tried searching this error to see if others may have faced the same error message, but unfortunately no luck. Any recommendations and or solutions would be greatly appreciated. The code and output can be found below.
kallisto quant -i Homo_sapiens.GRCh38.cdna.all -o test -t 8 --single -l 250 -s 30 SRR8668755_1M_subsample.fastq

[quant] fragment length distribution is truncated gaussian with mean = 250, sd = 30

Error: index input file could not be opened!

Usually errors like these occur when an input file doesn't exist or is corrupt. But you might also see these errors if you don't have permissions to access the file for whatever reason. Although unlikely, these can also occur if/when your system is out of memory.
Assuming the index is corrupt, I think just retry the command. For example:
curl -O ftp://ftp.ensembl.org/pub/release-108/fasta/homo_sapiens/cdna/Homo_sapiens.GRCh38.cdna.all.fa.gz

kallisto index \
    -i Homo_sapiens.GRCh38.cdna.all.release-108_k31.idx \
    Homo_sapiens.GRCh38.cdna.all.fa.gz

And test:
kallisto quant \
    -i Homo_sapiens.GRCh38.cdna.all.release-108_k31.idx \
    -o test \
    -t 8 \
    --single \
    -l 250 \
    -s 30 \
    SRR8668755_1M_subsample.fastq

To get this working using Windows/Anaconda Powershell, I just downloaded the cDNA sequences from https://ftp.ensembl.org/pub/release-108/fasta/homo_sapiens/cdna/ and moved them to my working directory. I also downloaded and extracted the Windows Kallisto executable v0.46.1 from https://pachterlab.github.io/kallisto/download.html. Then, with the FASTQ file in my working directory:
 .\kallisto_windows-v0.46.1\kallisto\kallisto.exe index `
    -i Homo_sapiens.GRCh38.cdna.all.release-108_k31.idx `
    .\Homo_sapiens.GRCh38.cdna.all.fa.gz

And tested with:
.\kallisto_windows-v0.46.1\kallisto\kallisto.exe quant `
    -i .\Homo_sapiens.GRCh38.cdna.all.release-108_k31.idx `
    -o test `
    -t 8 `
    --single `
    -l 250 `
    -s 30 `
    .\SRR8668755_1M_subsample.fastq


text:
Error installing VEP on mac
I have tried multiple different ways to get Ensemble's Variant Effect Predictor (VEP) to properly install on my mac, but I keep getting same error.
In /Users/username/ensembl-vep/htslib
make: Nothing to be done for `all'.
 - unpacking ./Bio/tmp/biodbhts.zip to ./Bio/tmp/
./Bio/tmp/Bio-DB-HTS-2.11 - moving files to ./biodbhts
 - making Bio::DB:HTS
clang-4.0: error: invalid version number in '-mmacosx-version-min=12.4'
Warning: ExtUtils::CBuilder not installed or no compiler detected
Proceeding with configuration, but compilation may fail during Build

Created MYMETA.yml and MYMETA.json
Creating new 'Build' script for 'Bio-DB-HTS' version '2.11'
Building Bio-DB-HTS
Error: no compiler detected to compile 'lib/Bio/DB/HTS.c'.  Aborting
ERROR: Shared Bio::DB:HTS library not found

--> When I check for the version of clang, I always get version 12, not sure where clang-4.0 is. I think error is specifically not correct path to current version of C++.

Do you have xcode utilities installed on your mac? What else have you tried?

You can try and get Perl's ExUtils::Builder, described here correctly installed and working, for example
perl -MCPAN -e shell
force install Bundle::CPAN

You might want to use brew to make the build safe.
I don't do this stuff, I use OSX, I would install ensemble's VEP via conda
conda install -c bioconda ensembl-vep

You might want to (good idea) create a specific environment and obviously download Anaconda or miniconda.
conda create --name vep
conda activate vep

It will work, the output is below. I didn't proceed with the installation because of the shear amount of stuff downloaded into a given conda environment.  There's everything in there but the kitchen sink, like there is mafft, t-coffee, clustalw, muscle and clustalo required in the download. I am honestly clueless because they are all alignment algorithms. Anyway I honestly wouldn't chase the MCPAN stuff, in this instance I definitely recommend conda. Finally, you honestly you don't want all this extra stuff on your hard drive when you've finished just delete the entire environment.
The final point I'd like to make is muscle5 (which is also in there) has recently significantly been modified making all other aligners obsolete on the runtime issue, but has to be run under the super5 mode.
## Package Plan ##

  environment location: /opt/anaconda3/envs/XXX # environment was edited  out

  added / updated specs:
    - ensembl-vep

The following packages will be downloaded:

    package                    |            build
    ---------------------------|-----------------
    argtable2-2.13             |    h1de35cc_1001         2.7 MB  conda-forge
    clustalo-1.2.4             |       h9722bc1_5         331 KB  bioconda
    clustalw-2.1               |       hcd10b59_7         326 KB  bioconda
    ensembl-vep-108.2          | pl5321h4a94de4_0         3.1 MB  bioconda
    expat-2.4.9                |       he9d5cce_0         127 KB
    htslib-1.14                |       hc38c3fb_0         2.0 MB  bioconda
    libdeflate-1.7             |       h9ed2024_5          52 KB
    libgfortran-5.0.0          |11_3_0_hecd8cb5_28         142 KB
    libgfortran5-11.3.0        |      h9dfd629_28         1.4 MB
    libzlib-1.2.12             |       hfd90126_3          65 KB  conda-forge
    mafft-7.508                |       ha5712d3_0         8.7 MB  bioconda
    muscle-5.1                 |       hb339e23_1         298 KB  bioconda
    mysql-connector-c-6.1.11   |       hccea1a4_1         1.2 MB
    paml-4.9                   |       hb4d813b_6         1.2 MB  bioconda
    perl-algorithm-diff-1.1903 | pl5321hdfd78af_3          36 KB  bioconda
    perl-base-2.23             | pl5321hdfd78af_2          12 KB  bioconda
    perl-bio-asn1-entrezgene-1.73| pl5321hdfd78af_3          24 KB  bioconda
    perl-bio-coordinate-1.007001| pl5321hdfd78af_3          36 KB  bioconda
    perl-bio-db-hts-3.01       | pl5321hc7384eb_4         116 KB  bioconda
    perl-bio-featureio-1.6.905 | pl5321hdfd78af_4          41 KB  bioconda
    perl-bio-samtools-1.43     | pl5321h1f540d2_3         161 KB  bioconda
    perl-bio-tools-phylo-paml-1.7.3| pl5321hdfd78af_3          58 KB  bioconda
    perl-bio-tools-run-alignment-clustalw-1.7.4| pl5321hdfd78af_3          25 KB  bioconda
    perl-bio-tools-run-alignment-tcoffee-1.7.4| pl5321hdfd78af_5          25 KB  bioconda
    perl-bioperl-1.7.8         |       hdfd78af_0           8 KB  bioconda
    perl-bioperl-core-1.7.8    | pl5321h527b516_0         1.8 MB  bioconda
    perl-bioperl-run-1.007003  | pl5321hdfd78af_0         340 KB  bioconda
    perl-business-isbn-3.007   | pl5321hdfd78af_0          18 KB  bioconda
    perl-business-isbn-data-20210112.006| pl5321hdfd78af_0          21 KB  bioconda
    perl-capture-tiny-0.48     | pl5321hdfd78af_2          19 KB  bioconda
    perl-class-data-inheritable-0.09| pl5321hdfd78af_0          10 KB  bioconda
    perl-constant-1.33         | pl5321hdfd78af_2          12 KB  bioconda
    perl-data-dumper-2.183     | pl5321ha5712d3_1          38 KB  bioconda
    perl-dbd-mysql-4.046       | pl5321hcd10b59_4          94 KB  bioconda
    perl-dbi-1.643             | pl5321ha5712d3_1         602 KB  bioconda
    perl-devel-stacktrace-2.04 | pl5321hdfd78af_1          17 KB  bioconda
    perl-digest-hmac-1.04      | pl5321hdfd78af_0          10 KB  bioconda
    perl-digest-md5-2.58       | pl5321ha5712d3_1          18 KB  bioconda
    perl-encode-locale-1.05    | pl5321hdfd78af_7          13 KB  bioconda
    perl-exception-class-1.45  | pl5321hdfd78af_0          21 KB  bioconda
    perl-file-listing-6.15     | pl5321hdfd78af_0          14 KB  bioconda
    perl-file-slurp-tiny-0.004 | pl5321hdfd78af_2          10 KB  bioconda
    perl-file-sort-1.01        | pl5321hdfd78af_3          20 KB  bioconda
    perl-file-spec-3.48_01     | pl5321hdfd78af_2           6 KB  bioconda
    perl-getopt-long-2.54      | pl5321hdfd78af_0          31 KB  bioconda
    perl-html-parser-3.80      | pl5321hb3ea9af_0          67 KB  bioconda
    perl-html-tagset-3.20      | pl5321hdfd78af_4          14 KB  bioconda
    perl-http-cookies-6.10     | pl5321hdfd78af_0          22 KB  bioconda
    perl-http-daemon-6.14      | pl5321hdfd78af_0          20 KB  bioconda
    perl-http-date-6.05        | pl5321hdfd78af_0          13 KB  bioconda
    perl-http-message-6.36     | pl5321hdfd78af_0          54 KB  bioconda
    perl-http-negotiate-6.01   | pl5321hdfd78af_4          14 KB  bioconda
    perl-io-html-1.004         | pl5321hdfd78af_0          16 KB  bioconda
    perl-io-socket-ssl-2.074   | pl5321hdfd78af_0         162 KB  bioconda
    perl-io-string-1.08        | pl5321hdfd78af_4          12 KB  bioconda
    perl-io-tty-1.16           | pl5321ha5712d3_1          30 KB  bioconda
    perl-ipc-run-20200505.0    | pl5321hdfd78af_0          76 KB  bioconda
    perl-libwww-perl-6.67      | pl5321hdfd78af_0          97 KB  bioconda
    perl-libxml-perl-0.08      | pl5321hdfd78af_3          40 KB  bioconda
    perl-lwp-mediatypes-6.04   | pl5321hdfd78af_1          24 KB  bioconda
    perl-mime-base64-3.16      | pl5321ha5712d3_2          16 KB  bioconda
    perl-net-http-6.22         | pl5321hdfd78af_0          21 KB  bioconda
    perl-net-ssleay-1.92       | pl5321h5d619eb_1         299 KB  bioconda
    perl-ntlm-1.09             | pl5321hdfd78af_5          16 KB  bioconda
    perl-perlio-gzip-0.20      | pl5321h1f540d2_3          16 KB  bioconda
    perl-sereal-4.019          | pl5321hdfd78af_0          11 KB  bioconda
    perl-sereal-decoder-4.025  | pl5321ha5712d3_1         156 KB  bioconda
    perl-sereal-encoder-4.025  | pl5321ha5712d3_1         365 KB  bioconda
    perl-set-intervaltree-0.12 | pl5321h9722bc1_2          31 KB  bioconda
    perl-socket-2.027          | pl5321ha5712d3_3          29 KB  bioconda
    perl-sub-uplevel-0.2800    | pl5321ha5712d3_4          14 KB  bioconda
    perl-test-deep-1.130       | pl5321hdfd78af_0          37 KB  bioconda
    perl-test-differences-0.69 | pl5321hdfd78af_0          19 KB  bioconda
    perl-test-exception-0.43   | pl5321hdfd78af_3          15 KB  bioconda
    perl-test-harness-3.44     | pl5321hdfd78af_0          85 KB  bioconda
    perl-test-most-0.38        | pl5321hdfd78af_0          20 KB  bioconda
    perl-test-warn-0.36        | pl5321hdfd78af_2          14 KB  bioconda
    perl-text-csv-2.01         | pl5321hdfd78af_0          70 KB  bioconda
    perl-text-diff-1.45        | pl5321hdfd78af_1          23 KB  bioconda
    perl-time-local-1.30       | pl5321hdfd78af_0          12 KB  bioconda
    perl-timedate-2.33         | pl5321hdfd78af_2          31 KB  bioconda
    perl-tree-dag_node-1.32    | pl5321hdfd78af_0          43 KB  bioconda
    perl-try-tiny-0.31         | pl5321hdfd78af_1          19 KB  bioconda
    perl-uri-5.12              | pl5321hdfd78af_0          59 KB  bioconda
    perl-url-encode-0.03       | pl5321h527b516_0          12 KB  bioconda
    perl-www-robotrules-6.02   | pl5321hdfd78af_4          14 KB  bioconda
    perl-xml-dom-1.46          | pl5321hdfd78af_1          59 KB  bioconda
    perl-xml-dom-xpath-0.14    | pl5321hdfd78af_2          12 KB  bioconda
    perl-xml-parser-2.44       | pl5321h50c80e2_8         162 KB  bioconda
    perl-xml-regexp-0.04       | pl5321hdfd78af_3          10 KB  bioconda
    perl-xml-xpathengine-0.14  | pl5321hdfd78af_3          28 KB  bioconda
    t-coffee-11.00.8cbe486     |       h26a2512_0         1.1 MB  bioconda
    unzip-6.0                  |       h3b3adb9_0          71 KB
    ------------------------------------------------------------
                                           Total:        28.7 MB

The following NEW packages will be INSTALLED:

  argtable2          conda-forge/osx-64::argtable2-2.13-h1de35cc_1001
  blast              bioconda/osx-64::blast-2.12.0-h0370960_3
  clustalo           bioconda/osx-64::clustalo-1.2.4-h9722bc1_5
  clustalw           bioconda/osx-64::clustalw-2.1-hcd10b59_7
  ensembl-vep        bioconda/noarch::ensembl-vep-108.2-pl5321h4a94de4_0
  entrez-direct      bioconda/osx-64::entrez-direct-16.2-h193322a_1
  expat              pkgs/main/osx-64::expat-2.4.9-he9d5cce_0
  htslib             bioconda/osx-64::htslib-1.14-hc38c3fb_0
  libdeflate         pkgs/main/osx-64::libdeflate-1.7-h9ed2024_5
  libgfortran        pkgs/main/osx-64::libgfortran-5.0.0-11_3_0_hecd8cb5_28
  libgfortran5       pkgs/main/osx-64::libgfortran5-11.3.0-h9dfd629_28
  libidn2            pkgs/main/osx-64::libidn2-2.3.2-h9ed2024_0
  libunistring       pkgs/main/osx-64::libunistring-0.9.10-h9ed2024_0
  libzlib            conda-forge/osx-64::libzlib-1.2.12-hfd90126_3
  llvm-openmp        pkgs/main/osx-64::llvm-openmp-14.0.6-h0dcd299_0
  mafft              bioconda/osx-64::mafft-7.508-ha5712d3_0
  muscle             bioconda/osx-64::muscle-5.1-hb339e23_1
  mysql-connector-c  pkgs/main/osx-64::mysql-connector-c-6.1.11-hccea1a4_1
  paml               bioconda/osx-64::paml-4.9-hb4d813b_6
  pcre               pkgs/main/osx-64::pcre-8.45-h23ab428_0
  perl               pkgs/main/osx-64::perl-5.32.1-0_h435f0c2_perl5
  perl-algorithm-di~ bioconda/noarch::perl-algorithm-diff-1.1903-pl5321hdfd78af_3
  perl-archive-tar   bioconda/noarch::perl-archive-tar-2.40-pl5321hdfd78af_0
  perl-base          bioconda/noarch::perl-base-2.23-pl5321hdfd78af_2
  perl-bio-asn1-ent~ bioconda/noarch::perl-bio-asn1-entrezgene-1.73-pl5321hdfd78af_3
  perl-bio-coordina~ bioconda/noarch::perl-bio-coordinate-1.007001-pl5321hdfd78af_3
  perl-bio-db-hts    bioconda/osx-64::perl-bio-db-hts-3.01-pl5321hc7384eb_4
  perl-bio-featureio bioconda/noarch::perl-bio-featureio-1.6.905-pl5321hdfd78af_4
  perl-bio-samtools  bioconda/osx-64::perl-bio-samtools-1.43-pl5321h1f540d2_3
  perl-bio-tools-ph~ bioconda/noarch::perl-bio-tools-phylo-paml-1.7.3-pl5321hdfd78af_3
  perl-bio-tools-ru~ bioconda/noarch::perl-bio-tools-run-alignment-clustalw-1.7.4-pl5321hdfd78af_3
  perl-bio-tools-ru~ bioconda/noarch::perl-bio-tools-run-alignment-tcoffee-1.7.4-pl5321hdfd78af_5
  perl-bioperl       bioconda/noarch::perl-bioperl-1.7.8-hdfd78af_0
  perl-bioperl-core  bioconda/osx-64::perl-bioperl-core-1.7.8-pl5321h527b516_0
  perl-bioperl-run   bioconda/noarch::perl-bioperl-run-1.007003-pl5321hdfd78af_0
  perl-business-isbn bioconda/noarch::perl-business-isbn-3.007-pl5321hdfd78af_0
  perl-business-isb~ bioconda/noarch::perl-business-isbn-data-20210112.006-pl5321hdfd78af_0
  perl-capture-tiny  bioconda/noarch::perl-capture-tiny-0.48-pl5321hdfd78af_2
  perl-carp          bioconda/noarch::perl-carp-1.38-pl5321hdfd78af_4
  perl-class-data-i~ bioconda/noarch::perl-class-data-inheritable-0.09-pl5321hdfd78af_0
  perl-common-sense  bioconda/noarch::perl-common-sense-3.75-pl5321hdfd78af_0
  perl-compress-raw~ bioconda/osx-64::perl-compress-raw-bzip2-2.201-pl5321h9722bc1_0
  perl-compress-raw~ bioconda/osx-64::perl-compress-raw-zlib-2.105-pl5321h9722bc1_0
  perl-constant      bioconda/noarch::perl-constant-1.33-pl5321hdfd78af_2
  perl-data-dumper   bioconda/osx-64::perl-data-dumper-2.183-pl5321ha5712d3_1
  perl-dbd-mysql     bioconda/osx-64::perl-dbd-mysql-4.046-pl5321hcd10b59_4
  perl-dbi           bioconda/osx-64::perl-dbi-1.643-pl5321ha5712d3_1
  perl-devel-stackt~ bioconda/noarch::perl-devel-stacktrace-2.04-pl5321hdfd78af_1
  perl-digest-hmac   bioconda/noarch::perl-digest-hmac-1.04-pl5321hdfd78af_0
  perl-digest-md5    bioconda/osx-64::perl-digest-md5-2.58-pl5321ha5712d3_1
  perl-encode        bioconda/osx-64::perl-encode-3.19-pl5321ha5712d3_1
  perl-encode-locale bioconda/noarch::perl-encode-locale-1.05-pl5321hdfd78af_7
  perl-exception-cl~ bioconda/noarch::perl-exception-class-1.45-pl5321hdfd78af_0
  perl-exporter      bioconda/noarch::perl-exporter-5.72-pl5321hdfd78af_2
  perl-exporter-tiny bioconda/noarch::perl-exporter-tiny-1.002002-pl5321hdfd78af_0
  perl-extutils-mak~ conda-forge/noarch::perl-extutils-makemaker-7.64-pl5321hd8ed1ab_0
  perl-file-listing  bioconda/noarch::perl-file-listing-6.15-pl5321hdfd78af_0
  perl-file-slurp-t~ bioconda/noarch::perl-file-slurp-tiny-0.004-pl5321hdfd78af_2
  perl-file-sort     bioconda/noarch::perl-file-sort-1.01-pl5321hdfd78af_3
  perl-file-spec     bioconda/noarch::perl-file-spec-3.48_01-pl5321hdfd78af_2
  perl-getopt-long   bioconda/noarch::perl-getopt-long-2.54-pl5321hdfd78af_0
  perl-html-parser   bioconda/osx-64::perl-html-parser-3.80-pl5321hb3ea9af_0
  perl-html-tagset   bioconda/noarch::perl-html-tagset-3.20-pl5321hdfd78af_4
  perl-http-cookies  bioconda/noarch::perl-http-cookies-6.10-pl5321hdfd78af_0
  perl-http-daemon   bioconda/noarch::perl-http-daemon-6.14-pl5321hdfd78af_0
  perl-http-date     bioconda/noarch::perl-http-date-6.05-pl5321hdfd78af_0
  perl-http-message  bioconda/noarch::perl-http-message-6.36-pl5321hdfd78af_0
  perl-http-negotia~ bioconda/noarch::perl-http-negotiate-6.01-pl5321hdfd78af_4
  perl-io-compress   bioconda/osx-64::perl-io-compress-2.201-pl5321h9722bc1_0
  perl-io-html       bioconda/noarch::perl-io-html-1.004-pl5321hdfd78af_0
  perl-io-socket-ssl bioconda/noarch::perl-io-socket-ssl-2.074-pl5321hdfd78af_0
  perl-io-string     bioconda/noarch::perl-io-string-1.08-pl5321hdfd78af_4
  perl-io-tty        bioconda/osx-64::perl-io-tty-1.16-pl5321ha5712d3_1
  perl-io-zlib       bioconda/noarch::perl-io-zlib-1.11-pl5321hdfd78af_0
  perl-ipc-run       bioconda/noarch::perl-ipc-run-20200505.0-pl5321hdfd78af_0
  perl-json          bioconda/noarch::perl-json-4.10-pl5321hdfd78af_0
  perl-json-xs       bioconda/osx-64::perl-json-xs-2.34-pl5321hcd10b59_5
  perl-libwww-perl   bioconda/noarch::perl-libwww-perl-6.67-pl5321hdfd78af_0
  perl-libxml-perl   bioconda/noarch::perl-libxml-perl-0.08-pl5321hdfd78af_3
  perl-list-moreuti~ bioconda/noarch::perl-list-moreutils-0.430-pl5321hdfd78af_0
  perl-list-moreuti~ bioconda/osx-64::perl-list-moreutils-xs-0.430-pl5321ha5712d3_1
  perl-lwp-mediatyp~ bioconda/noarch::perl-lwp-mediatypes-6.04-pl5321hdfd78af_1
  perl-mime-base64   bioconda/osx-64::perl-mime-base64-3.16-pl5321ha5712d3_2
  perl-net-http      bioconda/noarch::perl-net-http-6.22-pl5321hdfd78af_0
  perl-net-ssleay    bioconda/osx-64::perl-net-ssleay-1.92-pl5321h5d619eb_1
  perl-ntlm          bioconda/noarch::perl-ntlm-1.09-pl5321hdfd78af_5
  perl-parent        bioconda/noarch::perl-parent-0.236-pl5321hdfd78af_2
  perl-pathtools     bioconda/osx-64::perl-pathtools-3.75-pl5321ha5712d3_3
  perl-perlio-gzip   bioconda/osx-64::perl-perlio-gzip-0.20-pl5321h1f540d2_3
  perl-scalar-list-~ bioconda/osx-64::perl-scalar-list-utils-1.62-pl5321ha5712d3_0
  perl-sereal        bioconda/noarch::perl-sereal-4.019-pl5321hdfd78af_0
  perl-sereal-decod~ bioconda/osx-64::perl-sereal-decoder-4.025-pl5321ha5712d3_1
  perl-sereal-encod~ bioconda/osx-64::perl-sereal-encoder-4.025-pl5321ha5712d3_1
  perl-set-interval~ bioconda/osx-64::perl-set-intervaltree-0.12-pl5321h9722bc1_2
  perl-socket        bioconda/osx-64::perl-socket-2.027-pl5321ha5712d3_3
  perl-sub-uplevel   bioconda/osx-64::perl-sub-uplevel-0.2800-pl5321ha5712d3_4
  perl-test-deep     bioconda/noarch::perl-test-deep-1.130-pl5321hdfd78af_0
  perl-test-differe~ bioconda/noarch::perl-test-differences-0.69-pl5321hdfd78af_0
  perl-test-excepti~ bioconda/noarch::perl-test-exception-0.43-pl5321hdfd78af_3
  perl-test-harness  bioconda/noarch::perl-test-harness-3.44-pl5321hdfd78af_0
  perl-test-most     bioconda/noarch::perl-test-most-0.38-pl5321hdfd78af_0
  perl-test-warn     bioconda/noarch::perl-test-warn-0.36-pl5321hdfd78af_2
  perl-text-csv      bioconda/noarch::perl-text-csv-2.01-pl5321hdfd78af_0
  perl-text-diff     bioconda/noarch::perl-text-diff-1.45-pl5321hdfd78af_1
  perl-time-local    bioconda/noarch::perl-time-local-1.30-pl5321hdfd78af_0
  perl-timedate      bioconda/noarch::perl-timedate-2.33-pl5321hdfd78af_2
  perl-tree-dag_node bioconda/noarch::perl-tree-dag_node-1.32-pl5321hdfd78af_0
  perl-try-tiny      bioconda/noarch::perl-try-tiny-0.31-pl5321hdfd78af_1
  perl-types-serial~ bioconda/noarch::perl-types-serialiser-1.01-pl5321hdfd78af_0
  perl-uri           bioconda/noarch::perl-uri-5.12-pl5321hdfd78af_0
  perl-url-encode    bioconda/osx-64::perl-url-encode-0.03-pl5321h527b516_0
  perl-www-robotrul~ bioconda/noarch::perl-www-robotrules-6.02-pl5321hdfd78af_4
  perl-xml-dom       bioconda/noarch::perl-xml-dom-1.46-pl5321hdfd78af_1
  perl-xml-dom-xpath bioconda/noarch::perl-xml-dom-xpath-0.14-pl5321hdfd78af_2
  perl-xml-parser    bioconda/osx-64::perl-xml-parser-2.44-pl5321h50c80e2_8
  perl-xml-regexp    bioconda/noarch::perl-xml-regexp-0.04-pl5321hdfd78af_3
  perl-xml-xpatheng~ bioconda/noarch::perl-xml-xpathengine-0.14-pl5321hdfd78af_3
  t-coffee           bioconda/osx-64::t-coffee-11.00.8cbe486-h26a2512_0
  unzip              pkgs/main/osx-64::unzip-6.0-h3b3adb9_0
  wget               pkgs/main/osx-64::wget-1.21.3-h6dfd666_0

VEP requires gcc and g++ for installation1. Before running the build script, try setting the CC and CXX environment variables to the paths of your gcc and g++ compilers, respectively:
export CC=gcc
export CXX=g++

You'll need to install these compilers if they don't already exist on your system.


text:
split fast reference fast file by 22 chromosome
HG38 human reference is the collection of reference from all 22 chromosomes,x,y chromosomes, and some others. I was targeting to split the reference file into 22 chromosomes fasta files. I am following this link https://crashcourse.housegordon.org/split-fasta-files.html to split the fasta files.More specifically I want to split the fasta file into all possible (placed,unlocalized and unplaced) sequence and the merge placed and localized sequences.
However, I am unsure where to put the HLA part of the reference. Should I put it in Chromosome 6? Also, should I dump the unknown chromosomes as I am specifically interested in chromosomes 1-22?
Any suggestion on splitting the fasta file into 22 chromosomes would be nice.

"should I dump the unknown chromosomes as I am specifically interested in chromosomes 1-22". That's up to you. We cannot tell you what your data preservation strategy should be.

Is the question about the computing/parsing side (how do I make one file per sequence from a multi-fasta file?) or is it more about the biology of what the different chromosomes in hg38 actually are? You seem to be asking about both sides but I can't quite understand.

It is more about biology

Then please [edit] your question and make that clear. As it stands, I read this as asking about how to split a multi-fasta file. One thing though: if the people who generated the genome couldn't put the extra sequences in a specific chromosome, what makes you think you will be able to? This feels like an XY problem. Why do you think you need to end up with exactly 22 chromosomes? Can you explain what the final objective is here?

The HLA region is already a part of chromosome 6. When you split the file, it will just be included with the rest of the chromosome. Whether or not you include the remaining contigs depends on your specific objectives. Inclusion of the unplaced and unlocalized scaffolds usually helps to reduce mapping errors, for example. To split/extract only chr1-22 from the UCSC hg38.fa.gz, we can just use awk. With the following in a file called script.awk:
BEGIN {
    for(i=1;i<=22;i++) {
        arr["chr" i]
    }
}

/^>/ {
    f=substr($0,2)
}

f in arr {
    print > f ".fa"
}

Run using:
awk -f script.awk <(zcat hg38.fa.gz)

Of course, if that's what the OP wants, it would make more sense to just download the individual chromosomes from https://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/.


text:
Need to convert .bed to .vcf. Can the reference build (and necessary .fasta) be determined for a .bed file?
I have some .bed files (and .bim and .fam files) containing data for a number of different samples, and I need to convert them to .vcf.
I've found bed2vcf which is from the bedr R package, and also available as a stand-alone python program of the same name, but it requires a .fasta file corresponding to the reference genome used to create the .bed file. So we need to know what genome reference build was used to get the right .fasta.
Question Is there any tool or way to inspect the files I have and determine what genome reference build is required to do a correct conversion to .vcf, or is there any other method?
There is a bed2vcf related post here.

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

So it seems like you've (very reasonably) got confused between the two different types of bed file - one of them is used to store genomic coordinate data and then one is used as part of the plink data set. You want the latter, but I function you used is for the former. Annoying I know...

There are two different kinds of .bed files in bioinformatics; one related to storing genomic coordinate data and one which is part of the plink binary dataset, so it's important to check which kind of file you are actually converting. If you have .bed, .bim and .fam files, you could use PLINK 2.0 to export them as VCF:
plink2 --bfile [prefix] --export vcf --out <prefix>

The relevant options from the docs:
<...>

  --bfile  <prefix> ['vzs'] : Specify .bed + .bim[.zst] + .fam prefix

<...>

  --out <prefix>     : Specify prefix for output files.

<...>

  --export <output format(s)...> [{01 | 12}] ['bgz'] ['id-delim='<char>]
           ['id-paste='<column set descriptor>] ['include-alt']
           ['omit-nonmale-y'] ['spaces'] ['vcf-dosage='<field>] ['ref-first']
           ['bits='<#>] ['sample-v2']
    Create a new fileset with all filters applied.  The following output
    formats are supported:

...

    * 'vcf',     : VCF (default version 4.3).  If PAR1 and PAR2 are present,
      'vcf-4.2',   they are automatically merged with chrX, with proper
      'bcf',       handling of chromosome codes and male ploidy.
      'bcf-4.2'    When the 'bgz' modifier is present, the VCF file is
                   block-gzipped.  (This always happens with BCF output.)
                   The 'id-paste' modifier controls which .psam columns are
                   used to construct sample IDs (choices are maybefid, fid,
                   iid, maybesid, and sid; default is maybefid,iid,maybesid),
                   while the 'id-delim' modifier sets the character between the
                   ID pieces (default '_').

...

+1, but it's probably best to use plink2 for this, I think it's much quicker and can handle some different kind of variants I think.

Oh wow. I've never worked with PLINK so this is news to me. Why in the world would they call their format "bed"? That's just...

@user438383 Updated, thanks!

@terdon as an experienced bioinformatician this shouldn’t surprise you :)


text:
How can I resolve "missing name of population" message in PGDSpider?
I want to use Arlequin for AMOVA. I have my FASTA files which I want to convert into .arp files which are recognizable by arlequin. When I try to convert it using Population Genetics Data (format) Spider, or PGDSpider by Lisher and Excoffier (2018).
The error message says the name of my populations is missing.
Can anyone tell me what's wrong with my FASTA files? Check the following example:
>CcHr1
ATAAAGATATTGGTACGTTGTATATATTATTGGGTGTGTGGTGTGGTATAGTTGGTACAGGATTATCATT
ACTAATTCGTCTTGAGTTAGGTACTTCTGGTGTCTTGACTGATGAGCATTTTTATAATGTAATTGTAACA
GCCCACGCTTTTGTAATGATTTTTTTTATGGTTATACCGATTATAATTGGTGGTTTTGGAAATTGAATAG
TTCCTTTATTAATTGGTGCTCCAGACATAAGGTTCCCTCGGATAAATAATATAAGATTTTGATTATTACC
TCCTGCTTTTGTACTACTGATTAGTAGAAGCCTTGTTGAGGGAGGGGCTGGGACTGGCTGAACTGTTTAT
CCTCCACTTAGGAGGGTAACAGGACATAGTGGGGCTTCTGTAGACTTGGCTATTTTTTCCTTACATCTCG
CGGGCATATCATCAATTTTGGGGGCAATTAATTTTATTACAACGATCTTCAATATGCGGTCCCCTGGGTT
AACACTCGAGCGGATAAGGCTATTTGTATGGTCTATTTTAGTAACAGTTTTTTTATTGTTACTATCACTT
CCAGTTCTTGCAGGGGCAATTACGATGCTCTTAACTGATCGTAATTTTAGTACTTCATTTTTTGATCCAG
>CcMath1
ATAAAGATATTGGTACGTTATATATATTATTGGGTGTGTGGTGTGGAATAGTTGGCACAGGGTTGTCATT
ATTAATTCGTCTTGAGTTAGGTACTTCTGGTATCTTGACTGATGAGCATTTTTATAACGTAATTGTAACA
GCCCACGCTTTTGTAATAATTTTTTTTATGGTTATGCCAATTATAATTGGTGGTTTTGGAAATTGAATAG
TCCCTTTATTAATTGGTGCTCCAGACATAAGGTTCCCTCGAATAAATAATATAAGATTTTGACTATTACC
TCCTGCGTTTGTACTGCTGATTAGTAGAAGCCTTGTTGAGGGAGGGGCTGGAACTGGCTGAACTGTTTAT
CCTCCACTTAGAAGGGTAACAGGGCATAGTGGGGCTTCTGTAGACTTGGCGATCTTTTCCTTGCATCTCG
CGGGTATATCATCAATTTTGGGGGCAATTAATTTTATTACAACGATCTTCAATATGCGGTCCCCTGGGTT
AACACTCGAACGGATAAGGTTATTTGTATGGTCTATTTTAGTAACAGTTTTTTTATTGTTATTGTCACTT
CCAGTTCTTGCAGGGGCTATTACGATACTCTTAACTGATCGTAATTTTAGTACTTCGTTTTTTGACCCAG

As you can see I have to populations but it only names it as pop1.

Is this something I can do inside arlequin? Because so far I have only though of manipulating the flat file into distinguishing populations. This is a very tedious process I would like to avoid

You could try specifying the population keyword in the record header. From the PGDSpider manual (version 2.1.1.5 (May 2018)), at the bottom of page 67:

Population definition in header line: String after “population:” until
the first space is parsed as population name e.g.: >sample1 population:pop1 gene:HMA4

A good response.

It's important to understand what Arlequin is doing for MANOVA. It will compare two or more populations to assess whether there is a difference in the allele frequency between the populations. It will also assess "interaction" (interpreting that in "gene flow" needs some thinking).
The basic theory is that if two or more populations cannot be distinguished under a formal parametric test there must be very strong gene flow supporting the lack of population structure between these populations.
Geographic foci Population designation is usually done between geographic foci, thus population A is usually geographic location (focus) A and population B is a geographically separated region (focus). Thus the fasta header needs to know the geographic region of the sample. Thats it.
Population info The issue is that information is required within the fasta header. Thus population delineation is therefore essential. This is confirmed in the PGDSpider manual below the line for Arlequin: which you don't need to do because thats what the conversion does.
Arlequin therefore needs to group each member of the population into a (usually geographically seperated) population to perform MANOVA. It will require this information specifying on the fasta header line. At the moment the data is saying there is a single population, it thinks all the sequences are from the same geographic location.
Procedure Thus it's likely simply a case of assigning two or more broad geographic foci for your samples. Then going through the fasta file and designating each sequence into one of the geographic focus under the keyword "population".
Just modifying @Steve's reading of the manual
> CcHr1 population:geoFocusA gene:X
ATAAAGATATTGGTACGTTGTATATATTATTGGGTGTGTGGTGTGGTATAGTTGGTACAGGATTATCATT

>CcMath1 population:geoFocusB gene:X
 ATAAAGATATTGGTACGTTATATATATTATTGGGTGTGTGGTGTGGAATAGTTGGCACAGGGTTGTCATT

Then run PGDSpider
java -Xmx1024m -Xms512m -jar PGDSpider2-cli.jar -inputfile examples\example_Structure.fasta -inputformat FASTA -outputfile examples\output_Arlequin.arp -outputformat ARLEQUIN -spid examples\Structure_Arlequin.spid

You did attain the critical anomaly in the sequence format,

each line of a sequence should have fewer than 80 characters

Your rows have 70 residues.

Genetic structure (only required for AMOVA):

The genetic structure specifies the hierarchical genetic structure of the samples. It is possible to define groups of populations.
o start of the subsection: [[Structure]]
o name for the genetic structure (string within ” ”):  StructureName = “name”
o number of groups defined in the structure (int value): NbGroups = 5
o group definitions (list containing the names of the samples belonging to the group,
entered within braces “{ }”):
Id1 2 ACTCGGGTTCGCGCGC # the first pseudo-haplotype ACTCGGGCTCACGCGC # the second pseudo-haplotype

NbGroups=2
Group ={
  population1
}
Group ={
  population2
  population3
}

Advantage The advantage of MANOVA is you can have LOADS of populations and it will assess which have high levels of gene flow between them. Statistically it is fantastic for loads of different reasons (quite complicated reasons). The disadvantage of output of MANOVA is there will likely be  significant differences between many of the populations - possibly all of them. This is a result of the power of parametric statistics, but may not be particularly helpful for the analysis. If that is the outcome switching to Fst is the next part of the analysis (and a different question). Goodluck and you can always post back if you get stuck.