text:
How to analyse Orthofinder results
I read this paper (Colletotrichum shisoi sp. nov., an anthracnose pathogen of Perilla frutescens in Japan: molecular phylogenetic, morphological and genomic evidence). We ran Orthofinder and also created an Upset plot. Unfortunately, I do not understand how they calculated the following percentages in the "Conservation of proteins amongst Colletotrichum species" section:

A total of 11,914 orthogroups with two or more proteins were
identified (Fig. 4b). Of these, 7,950 groups (74.0% proteins from C.
shisoi, 63.0% from C. higginsianum, 78.3% from C. incanum, 74.2% from
C. tofieldiae and 73.7% from C. graminicola) were conserved in all
five species (Supplementary Tables S3 and S4)

How is it possible to calculate the values?
Thank you in advance

Ah this is a core-genome calculation. Very interesting study indeed. Agreed its not clear but I think @CorteZero has the right idea.

I think they simply calculate the percentage of proteins in each species that belongs to an orthogroup conserved in the five species.
There should be a tab separated text file called Orthogroups.GeneCount.tsv (at least in the Orthofinder version used in the paper). This file contains counts of the number of genes for each species in each orthogroup. You can use this file to get the total number of genes in each species associated with those 7950 orthogroups and then calculate its proportion compared to the total number of genes.
Hope this helps


text:
How can I use sequences with different lengths in arlequin?
I have mtDNA sequences of the cox1 gene from various populations of one species. I want to perform an AMOVA allele frequency analysis via Arlequin. This allele frequency package however requires my sequences to have the same length.

Is it possible to use my sequences?
Do I have to align them first maybe?

The length difference is most likely due to indels/deletions. Do I have to use specific parts of the sequence for each gene of the subunit?

Please [edit] your question and give us some more details. Are you talking about protein sequences? DNA? mRNA? Arlequin can perform various types of analysis, what analysis are you trying to do?

Thanks, that helps. I don't know what you mean by subunit though, do you mean subunit of the gene's protein product, of cytochrome c oxidase? Or are you working on some sort of protein complex formed by the protein products of multiple genes? In any case, yes, you would very likely need to align the sequences first. Why not at least try that and see if it works?

Oh BTW @Nickmofoe the correct name for the gene is cox1. Please avoid using the old nomenclature of COI, this can create confusion.

There are various workflows that do this that are easy to google, for example here. Can you clarify how these don't work for you?

All that is needed is to align the sequences within a formal aligner like MAFFT or muscle. Then pass that through the PDGSpider converter. Arlequin should just take it from there.
Manual annotation of the alignment is normal e.g. seqoTron and is traditionally regarded as the best approach rather than do this algorithmically. Thus I'm suggesting (strongly suggesting) do this prior the PGDSpider conversion. If an alignment has not been performed Arlequin should throw an error.
If there are trailing gaps (indels or '-') at the beginning and end of the alignment Arlequin will/should ignore these. You can of course simply cut them out with the alignment editor (Mega X will also work for alignment editing). The important point is that Arlequin must correctly identify homologues its an allele frequency package and nothing more. If there are errors in the alignment Arlequin will believe there is population structure when it could be absent.
I can't honestly believe that Arlequin will count an indel '-' as an informative site. It should just ignore those sites, because it is a very respectable package that has been around for years and years.
Okay, I am now signing off for Christmas... with occasional exceptions.

source: Stack Exchange

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

sort: votes


text:
How can I fix "Ill-defined genetic structure" Error message in Arlequin?
I am performing an  AMOVA analysis via Arlequin.
At this point, I have:

collected sequences of mtDNA from the cox1 gene,
aligned the different genes using clustalX
converted the files into .arp format.

I then altered the structure of my populations in Arlequin so that all the different populations are in one group. When I tried to run AMOVA I get the error message

"Ill-Defined genetic structure".

Any clues as to what I have to adjust?

If one group is one population this will not work with Arlequin AMOVA. There must be more than one population (I might have misread the question however).
There could several reasons for this output. The most likely is there are insufficient SNPs within the cox1 locus for MANOVA to identify a population structure with sufficient statistical robustness. What it requires are a minimum of 3 independent alleles (SNPs), really 3 genotypes. If the SNPs are linked on the same bit of DNA that doesn't count - thats just one genotype. What I suspect is happening is there is a dearth (not enough) genetic information within cox1 to perform MANOVA.
There could be other explanations, but this is the first issue to address. This isn't just Arlequin it is standard theory for any allele frequency based calculations and why this is important is more complicated to explain.

Please note, it is the run-up to Christmas .. okay put Christmas on hold. To respond to the comments ...
Firstly, if you hold the data set it is very easy to identify the number of genotypes. There are past posts on this herein. I am not inclined to describe alternative approaches to an algorithm I've still not yet published :-) ... however it is really easy and needs an empirical test. It is good to assess number of genotypes per se Arlequin will do this analysis ... finding that within the output I dunno.
Extending loci is always good and for mtDNA I would not recommend with 16S nor cox2. rRNA should be avoided in allele frequency calculations, its okay for taxonomic assignment (bacteria rather than mtDNA). cox2 is a small locus. I would recommend for mtDNA

cobB, formerly called cytB
intergenic region of mtDNA

cobB is a good target, it has a lot more selection on it for within species analysis than cox1. The intergenic region is rarely used, its a good locus because of the higher rates of mutation therein.

Would you suggest using a different genetic marker for this type of analyses? I do not have access to SNPs so the next best thing would be other mtDNA genes like COX2 and rRNA.

@Nickmofoe answered ... Merry Christmas BTW.

Alright so: The COX1 locus ended up not being the issue. I made some corrections to the flat files and they seemed to run with no issues. I got the results and everything seems to check out. I concluded that the issue was not the locus itself but rather the incredibly complicated .arp format. Thanks for the advice though! Much appreciated and merry Christmas!!!

Merry Christmas to you @nickmofoe! The  "complicated format" was why Arlequin was never popular, because everybody felt it was unnecessarily complicated. This why your original question on PGDSpider was interesting because it as overdue requirement to make Arlequin accessible.

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


text:
DESeq2 for large number of samples takes too much RAM
I am trying to run a very large number of transposase-accessible chromatin (ATAC)-seq samples through DESeq2 to find peaks of differential chromatin accessible across my study genome.
I have about 100 samples, and 1,000,000 "genes" which are really peaks.
Doing this takes a significant amount of RAM; right now I am running on a 1 TB bigmem node and still getting an out-of-memory error. The code works perfectly fine for a small subset of samples.
Does anyone have an alternative strategy to use DESeq2 on a very large dataset? Should I look into doing a Wilcox-rank test myself? Thanks.

I'll preface this by saying that I don't think DESeq2 is the right tool to use for ATAC-Seq data. My own study of ATAC-Seq patterns [admittedly only a couple of runs that were our first exploration of the technique] indicates that it has a fractal expression which is difficult to properly summarise at any scale, especially at a bulk sample level. I found this poster from a quick search on Twitter, which recommends AIAP for ATAC-Seq analysis (but don't read too much into that; I have no experience with doing ATAC-Seq analysis). What I do have a lot of experience with is DESeq2.
DESeq2 has certain assumptions about the expression patterns of genes which are built into its differential expression model, and the further you get from gene-level analysis, the less appropriate DESeq2 is for data analysis (for example, Swish is a better approach for isoform level and allelic analysis).
However, if you're absolutely set on using DESeq2...  the DESeq2 data analysis documentation has a few different suggestions for speeding up analysis:

The above steps should take less than 30 seconds for most analyses. For experiments with complex designs and many samples (e.g. dozens of coefficients, ~100s of samples), one may want to have faster computation than provided by the default run of DESeq. We have two recommendations:

By using the argument fitType="glmGamPoi", one can leverage the faster NB GLM engine written by Constantin Ahlmann-Eltze. Note that glmGamPoi’s interface in DESeq2 requires use of test="LRT" and specification of a reduced design.
One can take advantage of parallelized computation. Parallelizing DESeq, results, and lfcShrink can be easily accomplished by loading the BiocParallel package, and then setting the following arguments: parallel=TRUE and BPPARAM=MulticoreParam(4), for example, splitting the job over 4 cores. However, some words of advice on parallelization: first, it is recommend to filter genes where all samples have low counts, to avoid sending data unnecessarily to child processes, when those genes have low power and will be independently filtered anyway; secondly, there is often diminishing returns for adding more cores due to overhead of sending data to child processes, therefore I recommend first starting with small number of additional cores. Note that obtaining results for coefficients or contrasts listed in resultsNames(dds) is fast and will not need parallelization. As an alternative to BPPARAM, one can register cores at the beginning of an analysis, and then just specify parallel=TRUE to the functions when called.

Also, applying a low count filter can help a lot.

Pre-filtering
While it is not necessary to pre-filter low count genes before running the DESeq2 functions, there are two reasons which make pre-filtering useful: by removing rows in which there are very few reads, we reduce the memory size of the dds data object, and we increase the speed of the transformation and testing functions within DESeq2. It can also improve visualizations, as features with no information for differential expression are not plotted.
Here we perform a minimal pre-filtering to keep only rows that have at least 10 reads total. Note that more strict filtering to increase power is automatically applied via independent filtering on the mean of normalized counts within the results function.

keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]

Alternatively, a popular filter is to ensure at least X samples with a count of 10 or more, where X can be chosen as the sample size of the smallest group of samples:

# unevaluated chunk...
keep <- rowSums(counts(dds) >= 10) >= X
dds <- dds[keep,]

Awesome, thanks for your input! Pre-filtering peaks that have a minimum count significantly reduces the space complexity. Of course, the million dollar question is: what minimum number of counts should each peak have... 5? 10? 100? Anyway, I think that's a question I should try and answer on my own. Thank you very much!

Adding to the previous comment, yes I am using parallel compute. I was not, however, using the NB GLM model; I am trying it right now and will report the speedup. w.r.t. using DESeq2 for ATAC-seq samples... I totally agree that the implicit assumptions make it a far-from-perfect algorithm to model the distribution in question. The only reason to use it is that other publications use it, which is a terrible reason. I think the Wilcoxon rank-sum test might be more appropriate here.

Absolutely, look at Wilcox-rank, not purely because of the computational overhead, but the parametric assumption that DESeq2 uses  are not correct on large data sets. This might describe your data sets too:

Abstract
DESeq2 and edgeR, have unexpectedly high false discovery rates.
Expanding the analysis to limma-voom, NOISeq, dearseq, and Wilcoxon rank-sum test, we found that [false discovery rates] FDR control is often failed except for the Wilcoxon rank-sum test
Particularly, the actual FDRs of DESeq2 and edgeR sometimes exceed 20% when the target FDR is 5%. Based on these results, for population-level RNA-seq studies with large sample sizes, we recommend the Wilcoxon rank-sum test.
Results and Discussion
DESeq2 and edgeR had only an 8% overlap in the DEGs they identified

I feel it's unfair to write that without mentioning the DESeq author response: in short, there are technical confounders that need to be controlled for. Wilcoxon gets fewer False Positives because of its lower power (while DESeq and edgeR are designed to maximize power), but the right thing to do is really to correct for confounding, whatever test you do next. That being said for this question Wilcoxon is probably reasonable (though a correction with RUV/SVA/... would be appropriate).

It is not a matter of parametric assumptions. The linked paper, as said already did not address the fact that in larger cohorts you often find samples that form clusters by some unwanted confounder nested with the variable you're testing. That is not unusual especially in human samples where you have all sorts of confounders, be it eating habits, drug consumption, dietary status etc that one often does not know of. Once this was addressed (RUVseq by the DESeq2 author) the main claim of that paper collapsed. In any case, at larger n I would use linear approaches such as limma-voom.

Hi @ATpoint thanks its a good discussion, albeit maybe not this site.

@M__, I think these comments are necessary to put answers into context. The paper you reference is basically an advertisement of the non-parametric tool the authors created. It's in the nature of these sort of papers that advertised tool is best and all others are not as good. That goes for any benchmark paper written by method developers, be it aligners, variant callers or stats frameworks. I am not judging this, it is just an observation. Like any of these papers, results should be read with a grain of salt and I think it is important to notice that this was not an unbiased benchmark.

@ATpoint we can have this discussion via the site  that you are famous on. There been an objection raised on this site and we must be mindful of that. I will raise this as a question on your site not quite now because it's Christmas Day (quite literally), but shortly. I agree 100% is an important topic and I'll explain my interest in the subject (I don't really do eukaryotes BTW).


text:
Read Clustal file in Python
This question was also asked on StackOverflow
I have a multiple sequence alignment (MSA) file derived from mafft in clustal format which I want to import into Python and save into a PDF file. I need to import the file and then highlight some specific words. I've tried to simply import the pdf of the MSA but after the highlight command doesn't work.
I need to print the file like this:

Multi.txt :
CLUSTAL format alignment by MAFFT FFT-NS-i (v7.453)

Consensus       --------------------------------------acgttttcgatatttatgccat
AMP             tttatattttctcctttttatgatggaacaagtctgcgacgttttcgatatttatgccat
                                                      **********************

Consensus       atgtgcatgttgtaaggttgaaagcaaaaatgaggggaaaaaaaatgaggtttttaataa
AMP             atgtgcatgttgtaaggttgaaagcaaaaatgaggggaaaaaaaatgaggtttttaataa
                ************************************************************

Consensus       ctacacatttagaggtctaggaaataaaggagtattaccatggaaatgtatttccctaga
AMP             ctacacatttagaggtctaggaaataaaggagtattaccatggaaatgtaattccctaga
                ************************************************** *********

Consensus       tatgaaatattttcgtgcagttacaacatatgtgaatgaatcaaaatatgaaaaattgaa
AMP             tatgaaatattttcgtgcagttacaacatatgtgaatgaatcaaaatatgaaaaattgaa
                ************************************************************

Consensus       atataagagatgtaaatatttaaacaaagaaactgtggataatgtaaatgatatgcctaa
AMP             atataagagatgtaaatatttaaacaaagaaactgtggataatgtaaatgatatgcctaa
                ************************************************************

Consensus       ttctaaaaaattacaaaatgttgtagttatgggaagaacaaactgggaaagcattccaaa
AMP             ttctaaaaaattacaaaatgttgtagttatgggaagaacaaactgggaaagcattccaaa
                ************************************************************

Consensus       aaaatttaaacctttaagcaataggataaatgttatattgtctagaaccttaaaaaaaga
AMP             aaaatttaaacctttaagcaataggataaatgttatattgtctagaaccttaaaaaaaga
                ************************************************************

Consensus       agattttgatgaagatgtttatatcattaacaaagttgaagatctaatagttttacttgg
AMP             agattttgatgaagatgtttatatcattaacaaagttgaagatctaatagttttacttgg
                ************************************************************

Consensus       gaaattaaattactataaatgttttattataggaggttccgttgtttatcaagaattttt
AMP             gaaattaaattactataaatgttttattataggaggttccgttgtttatcaagaattttt
                ************************************************************

Consensus       agaaaagaaattaataaaaaaaatatattttactagaataaatagtacatatgaatgtga
AMP             agaaaagaaattaataaaaaaaatatattttactagaataaatagtacatatgaatgtga
                ************************************************************

Consensus       tgtattttttccagaaataaatgaaaatgagtatcaaattatttctgttagcgatgtata
AMP             tgtattttttccagaaataaatgaaaatgagtatcaaattatttctgttagcgatgtata
                ************************************************************

Consensus       tactagtaacaatacaacattgga----------------------------------
AMP             tactagtaacaatacaacattggattttatcatttataagaaaacgaataataaaatg
                ************************

How can I import the alignment and print in the new PDF with the right alignment of the sequences.
Thanks

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.

What words do you need to highlight?

Ok, figured out a way, not sure its the best one,
nedd to install fpdf2 (pip install fpdf2)
from io import StringIO

from Bio import AlignIO  # Biopython 1.80

from fpdf import FPDF # pip install fpdf2

alignment = AlignIO.read("Multi.txt", "clustal")

stri = StringIO()

AlignIO.write(alignment, stri, 'clustal' )

# print(stri.getvalue())

stri_lines = [ i for i in stri.getvalue().split('\n')]

# print(stri_lines)

pdf = FPDF(orientation="P", unit="mm", format="A4")

# Add a page
pdf.add_page()

pdf.add_font('FreeMono', '', 'FreeMono.ttf')

pdf.set_font("FreeMono", size = 8)

for x in stri_lines:
    pdf.cell(0, 5, txt = x, border = 0,  new_x="LMARGIN" , new_y="NEXT",  align = 'L', fill = False)

    # print(len(x))

pdf.output("out.pdf")

output pdf  out.pdf :

Not sure why the file Header is changed, think is something within Biopythion (!!! ???) you can check adding:
with open('file_output.txt', 'w') as filehandler:
    AlignIO.write(alignment, filehandler, 'clustal')

I had to place a tff FreeMono font (Mono spaced font in my script directory) see: pdf.add_font('FreeMono', '', 'FreeMono.ttf') otherwise the alignement won't be printed in the correct way
Which fonts have the same width for every character?.
Attached a png of my pdf. See that you can highlight it,
using:
pdf.set_fill_color(255, 255, 0)

filling = False
for x in stri_lines:
    if 'Consensus' in x:
        filling = True
    else:
        filling = False
    pdf.cell(0, 5, txt = x, border = 0,  new_x="LMARGIN" , new_y="NEXT",  align = 'L', fill = filling)

or something similar you can highlight while printing:

pyMSAviz may help with this problem.
pyMSAviz is a Python tool with a convenient CLI/API implementation for visualizing MSA files.
For example, you can output MSA figure in PDF format with the following command.
pymsaviz -i Multi.txt -o msa_result.pdf --format clustal \
         --wrap_length 60 --show_count --color_scheme Identity

You can also use pyMSAviz's API to configure the appearance of the figure in more detail, such as highlighting.

Thanks for your solution Moshi. I personally think @pippo1980's solution is cooler because they leverage AlignIO within the StringIO (or StringIO within AlignIO). This permits a lot of flexibility in possible solutions because we use BioPython (AlignIO) all the time and is very cool OO.

personally I think your tool is supercool, spent last 3 days googling for such thing, but wasnt able to find anything ( I mean espript3 is supersupercool like jalview but no python API ... could be wrong). Wasted a day trying to understand how to get PyMuPDF  to find words, like 'ggaag' in my out.pdf with no luck. I am intrigued by your tool. Is there a way to color conserved residues all blue and non conserved like red ?? Was thinking about how to realize such a thing, Pyfpdf2 can highligh single cell and I should find a way to print each single alignment letter in a single cell and devising

a strategy to do that. Sorry I cannot upvote more than once.

yep found  --color_scheme Identity in your examples !!!

Hi @pippo1980. I'm pretty sure you could get yours working. Note that MEGA X alignment editor will colour conserved residues (but you don't get to choose the colour) and can't do non-conserved residues. Merry Christmas!

adding, a second answer because the OP request about

I need to import the file and then highlight some specific word

was making me uneasy. I kept using the Pyfpdf2 library and made use of the Biopython SeqRecord.letter_annotations attribute to track Nucleotides index composing the words found by RegEx. Tried to use the page.search_for() of the PyMuPDF library directly on the output pdf of my previous answer but was not getting word accross line breaks [if anybody know better, please help].
Here my code, it's a two half problem, first find the word (in my example is 'atat'), then write everything to pdf, I stole it from from /master/Bio/AlignIO/ClustalIO.py/   def write_alignment(self, alignment): from class ClustalWriter(SequentialAlignmentWriter) with some modding [I believe but not absolutely sure that ...A cell by definition is uniform in font and color... ]:
import re

from Bio import AlignIO  # Biopython 1.80

from fpdf import FPDF # pip install fpdf2

alignment = AlignIO.read("Multi.txt", "clustal")

word_to_find = 'atat'

for rec in alignment:

    rec.letter_annotations['highlight'] = 'N'*len(rec.seq)

# for rec in alignment:
#     print('rec.id : ', rec.id, 'highlight : ', rec.letter_annotations['highlight'])

for rec in alignment:

    # print(str(rec.seq))

    found = [m.start() for m in re.finditer( '(?='+re.escape(word_to_find)+')' , re.escape(str(rec.seq.replace('-', 'Z'))))]

    # print(rec.id , ' found : ', found)

    for i in found :

            a = list(rec.letter_annotations['highlight'])

            a[i:i+len(word_to_find)] = 'Y'*len(word_to_find)

            rec.letter_annotations['highlight'] = ''.join(i for i in a)

            # print(rec.letter_annotations['highlight'])

max_length = alignment.get_alignment_length()

# stolen from /master/Bio/AlignIO/ClustalIO.py/   def write_alignment(self, alignment): from class ClustalWriter(SequentialAlignmentWriter):
def write_alignment(alignment):

        x = 2

        y = 2

        pdf = FPDF(orientation="P", unit="mm", format="A4")

        # Add a page
        pdf.add_page()

        pdf.add_font('FreeMono', '', 'FreeMono.ttf')

        pdf.set_font("FreeMono", size = 12)

        pdf.set_fill_color(255, 255, 0)

        """Use this to write (another) single alignment to an open file."""
        if len(alignment) == 0:
            raise ValueError("Must have at least one sequence")
        if alignment.get_alignment_length() == 0:
            # This doubles as a check for an alignment object
            raise ValueError("Non-empty sequences are required")

        # Old versions of the parser in Bio.Clustalw used a ._version property
        try:
            version = str(alignment._version)
        except AttributeError:
            version = ""
        if not version:
            version = "1.81"
        if version.startswith("2."):
            # e.g. 2.0.x
            output = f"CLUSTAL {version} multiple sequence alignment"
        else:
            # e.g. 1.81 or 1.83
            output = f"CLUSTAL X ({version}) multiple sequence alignment"

        pdf.cell(0, y, txt = output , border = 0,  new_x='LMARGIN'  , new_y="NEXT" ,  align = 'L', fill = False)
        pdf.cell(0, 4*y, txt = '' , border = 0,  new_x='LMARGIN'  , new_y="NEXT" ,  align = 'L', fill = False)

        cur_char = 0
        # # max_length = len(alignment[0])
        max_length = alignment.get_alignment_length()

        if max_length <= 0:
            raise ValueError("Non-empty sequences are required")

        if "clustal_consensus" in alignment.column_annotations:
            star_info = alignment.column_annotations["clustal_consensus"]
        else:
            try:
                # This was originally stored by Bio.Clustalw as ._star_info
                star_info = alignment._star_info
            except AttributeError:
                star_info = None

        # keep displaying sequences until we reach the end

        rec_cnt = 0

        while cur_char != max_length:
            # calculate the number of sequences to show, which will
            # be less if we are at the end of the sequence
            if (cur_char + 50) > max_length:
                show_num = max_length - cur_char
            else:
                show_num = 50

            # go through all of the records and print out the sequences
            # when we output, we do a nice 80 column output, although this
            # may result in truncation of the ids.
            for record in alignment:

                rec_cnt += 1

                # Make sure we don't get any spaces in the record
                # identifier when output in the file by replacing
                # them with underscores:
                line = record.id[0:30].replace(" ", "_").ljust(50)

                pdf.cell(50, y, txt = line, border = 0,  new_x='RIGHT' , align = 'L', fill = False)

                line_seq = list(str(record.seq[cur_char : (cur_char + show_num)]))

                line_seq_highlight = list(str(record.letter_annotations['highlight'][cur_char : (cur_char + show_num)]))

                for i in range(len(line_seq)):

                    if line_seq_highlight[i] == 'N' :

                        filling = False

                    if line_seq_highlight[i] == 'Y' :

                        filling = True

                    else:

                        filling = False

                    pdf.cell(x, y, txt = line_seq[i], border = 0,  new_x= 'RIGHT' ,  align = 'C', fill = filling)

                pdf.cell(x, 2*y, txt = ' ' , border = 0,  new_x="LMARGIN" , new_y="NEXT" ,  align = 'L', fill = False)

                # now we need to print out the star info, if we've got it

                if star_info and rec_cnt == len(alignment):

                    pdf.cell(x, 2*y, txt = ' ' , border = 0,  new_x="LMARGIN" , new_y="NEXT" ,  align = 'L', fill = False)

                    rec_cnt = 0

                    # print(star_info)

                    line = (' ' * 50)

                    pdf.cell(50, y, txt = line, border = 0,  new_x='RIGHT' , align = 'L', fill = False)

                    star_info_seq = list(str(star_info[cur_char : (cur_char + show_num)]))

                    for i in range(len(star_info_seq)):

                        pdf.cell(x, y, txt = star_info_seq[i], border = 0,  new_x= 'RIGHT' ,  align = 'C', fill = False)

                pdf.cell(x, 2*y, txt = ' ' , border = 0,  new_x="LMARGIN" , new_y="NEXT" ,  align = 'L', fill = False)

            cur_char += show_num

            pdf.cell(x, 3*y, txt = ' ' , border = 0,  new_x="LMARGIN" , new_y="NEXT" ,  align = 'L', fill = False)

        pdf.output("out.pdf")

write_alignment(alignment)

I tested it, with the Multi.txt Clustal alignment from above, output is
out.pdf :

.
One last note:
needed to replace - with Z in the Consensus string to be able to get reliable results from the RegEX line : found = [m.start() for m in re.finditer( '(?='+re.escape(word_to_find)+')' , re.escape(str(rec.seq.replace('-', 'Z'))))] , not sure what is going on here, got the snippet from a SO post.


text:
Expending UpsetR plot with Orthofinder families and gene numbers
I can create an Upset plot with the input file, which can be found here
> library("UpSetR")
> orthogroups_df<- read.table("Orthogroups.GeneCount.tsv.txt",  header=T, stringsAsFactors = F)
> #All species
> selected_species <- colnames(orthogroups_df)[2:(ncol(orthogroups_df) -1)]
> selected_species
 [1] "Atha" "Vvin" "Slyc" "Stub" "Cann" "Natt" "Ntab" "Ntom" "Nsyl" "Ngla" "NQLD" "Nlab"
> head(orthogroups_df)
  Orthogroup Atha Vvin Slyc Stub Cann Natt Ntab Ntom Nsyl Ngla NQLD Nlab Total
1  OG0000000    0    0    0    0    0    0    0    0    0    0  965    3   968
2  OG0000001    0    0    0    0    1    0    0    0    0    0    3  448   452
3  OG0000002    0    0    0    0    1    0    0    1    1    0  313  120   436
4  OG0000003    0   26   25   39   93   21   63   36   33   46   15   16   413
5  OG0000004    1    0    9    4   42   34  154   11    8  109    2    6   380
6  OG0000005    0    0   20    1    2    1   70   43   91   34   61   44   367
> ncol(orthogroups_df)
[1] 14
> orthogroups_df[orthogroups_df > 0] <- 1
> upset(orthogroups_df, nsets = ncol(orthogroups_df), sets = rev(selected_species), keep.order = T, order.by = "freq")

I read this paper, and in Figure 2b and Supplementary Table S4, they had the gene numbers additionally to a number of families.

How can the upset plot be reproduced with gene numbers or/and create such a table?
Thank you in advance,


text:
bsub with -M and -hl a good idea?
I have been submitting my snakemake pipeline to a node-shared high performance computer (HPC) via bsub the following way:
bsub -W {resources.runtime_min} -n {threads} -M {resources.mem_mb} -R 'rusage[mem={resources.mem_mb}]' -o run_logs/{name}_%J.out -e run_logs/{name}_%J.err"

Works great, no problem.
However, I noticed that the memory limit was tasks*memory given. I am not aware of any memory issues on the HPC. The "issue" was not noticed until I asked why my memory limit was so high compared to what I gave with the bsub flag -M (I wanted to understand how everything worked).
I was recommend to use a -hl flag to hard cap the usage. It was actually an administator that recommended it, after I asked him about how it worked. They now want to implement it for other users. I found this in the documentation too, so I understand this.
But I see no one else is using this and it is also not suggested when you want to use the bsub outfiles to understand how much resources you need (how much memory should I use).
Does anyone have thoughts on this?
Somehow, it makes sense to use -hl, but I am hesitant since no one seems to use -hl for some reason, and I feel like I do not fully understand it.
I'm asking here, because I felt I kept bugging the administrator with the same questions :'-). Thank you in advance for any advice.

What this question appears to be about is capping the RAM usage on a node within a shared HPC.
Being 'nice' on an HPC is very honourable indeed, it is however the system administrators jurisdiction. It simply requires an email to them and they will provide the authoritative answer. Thus if there were RAM usage issues within the nodes then they will be aware of it. Limiting RAM on a job is risky because the biggest 'fall-over' in bioinformatics is a job running into a RAM bottleneck. If that happens the process just freezes and literally will take 1 log more run-time. If the machine doesn't have enough RAM for the job, thats a separate issue and easily solvable via cloud computing.
The individual who made the suggestion might be running into RAM bottleneck issues on the machine. However, they should inform the HPC administrator who could then provide an assessment whether all other users should cap their RAM usage. There are alternatives - the user experiencing RAM bottlenecks can parallelise their job and this will overcome RAM issues. This can be a bit tricky on a course-grain parallelisation with a shared HPC.
So in summary, I wouldn't do this personally and if I was concerned I would seek the advice of the administrator who will give the definitive response. If admin says 'do it' then do it, but it your jobs take way-way longer than they should I'd relax the constraint and inform admin.
Admin could enforce this constraint without the user base being aware of what was happening if it was a major issue on the server.
These days if there are HPC issues its easier to shift on a cloud computer and is not worth the arguments when everyone is trying to get their jobs done before the vacation.
Merry Christmas!

From the comments ...

System admin recommended this ... okay they should now inform the rest of the user base (which is what you are implying), unless there are been a specific issue. Thats a bit different and it is polite to comply with their request;
What you can do is justify to the system admin why you job requires extra RAM;
Systems admin might have automatically restricted the RAM access on the account - you can ask them to lift that;

You said,

"1 log more run-time"?

That is extremely important to understand. Bioinformatic jobs can use loads of RAM, if the job exceeds the RAM capacity of the node essentially the job "freezes". What it really does is slow the time to completion of the job dramatically, so I'm saying in my experience it's about 10-fold longer to complete than normal. That is what I am calling a RAM-bottleneck.
Final point
On the bsub parallelisation appears to be efficiently implemented and if RAM issues are occurring parallelising the job is the best solution. This would be a work around if the systems admin wants to restrict RAM usage and as a user you are wanting quick run-time (fast completion). Not all algorithms will parallelise but most will these days.

Thank you for the great answer @M__ ! Merry x-mas!


text:
How does one account for both single end and paired end reads as input in a Nextflow pipeline?
Suppose you have a pipeline that takes reads as input from the user and trims them.
// Declare syntax version
nextflow.enable.dsl=2

/*
 * pipeline input parameters
 */

params.reads = "$projectDir/data/*_{1,2}.fastq"

/*
 * 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)

}

This particular one takes paired end reads as input. What if the user wanted to trim both single end and paired end reads, how to account for single end reads as the input?
Best wishes.

Oh yes, I do upvote every response. and after implementing I accept one of the solutions.

Thanks @pubserfted and Merry Christmas!

@M__ Merry Christmas! Hope you have a good one. ^-^

Thank you @pubsurfted! I'm worked to the bone and will definitely enjoy the vacation.

Having to generate a sample sheet is just wasting time when you can detect it from the directory contents.

There are many ways to accomplish this. Assuming you don't need to process single-end and paired-end reads simultaneously (i.e. within the same run), one way using the fromFilePairs factory method is to just use an additional flag to define the number of files each emitted item is expected to hold. And then adjust your process definition accordingly:
params.reads = "${projectDir}/data/*_{1,2}.fastq.gz"
params.singleEnd = false

process TRIM {

    tag "${pair_id}"

    cpus 1

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

    output:
    tuple val(pair_id), path("${pair_id}/*.atria.fastq{,.gz}")

    script:
    def single = reads instanceof Path

    def read1 = !single ? /-r "${reads[0]}"/ : /-r "${reads}"/
    def read2 = !single ? /-R "${reads[1]}"/ : ''

    """
    atria \\
        -t ${task.cpus} \\
        ${read1} \\
        ${read2} \\
        -l 12 \\
        -q 30 \\
        -o "${pair_id}"
    """
}

workflow {

    Channel
        .fromFilePairs( params.reads, size: params.singleEnd ? 1 : 2 )
        .ifEmpty { error "Could not find any reads matching pattern: ${params.reads}" }
        .set { reads }

    TRIM( reads ).view()
}

Although we've used an additional flag, it does let the pipeline fail immediately if it can't find the files it expects. I think this approach is more explicit and provides a better user experience:
$ nextflow run main.nf --reads './data/*.fastq'
N E X T F L O W  ~  version 22.10.0
Launching `main.nf` [nasty_allen] DSL2 - revision: 2cbc311974
Could not find any reads matching pattern: ./data/*.fastq

 -- Check script 'main.nf' at line: 38 or see '.nextflow.log' file for more details

Doh! I forgot that these are single ended reads:
$ nextflow run main.nf --reads './data/*.fastq' --singleEnd
N E X T F L O W  ~  version 22.10.0
Launching `main.nf` [clever_gutenberg] DSL2 - revision: 2cbc311974
executor >  local (3)
[8e/f90c90] process > TRIM (X) [100%] 3 of 3 ✔
[Y, /path/to/work/e6/48df1485bb7de8355ed8ce0f1b53da/Y/Y.atria.fastq]
[Z, /path/to/work/eb/aca3389344ac03ff283ded0c46a7a6/Z/Z.atria.fastq]
[X, /path/to/work/8e/f90c90210103c03a3fbde697e01227/X/X.atria.fastq]

My personal strategy is to always read data from a samplesheet which then auto-detects whether a sample is paired-end or single-end and stores this information in a so-called meta map (inspired by nf-core) which then passes this information to the modules.
Example:
Say you have a samplesheet that is a three-column csv file with a header:
$ cat samplesheet.csv
sample,r1,r2
sample1,/Users/atpoint/sample1_R1.fq.gz,/Users/atpoint/sample1_R2.fq.gz
sample2,/Users/atpoint/sample2_R1.fq.gz,

Read this into Nextflow and detect whether r2 has an entry. If so it is paired-end, else it is single-end:

# main.nf

nextflow.enable.dsl = 2

// Channel for the samplesheet
ch_samplesheet = Channel.fromPath("/Users/atpoint/samplesheet.csv")

// Parse it line by line
ch_reads = ch_samplesheet.splitCsv(header:true).map {

    // This is the read1 and read2 entry
    r1 = it['r1']
    r2 = it['r2']

    // Detect wiether single-end or paired-end
    is_singleEnd = r2.toString()=='' ? true : false

    // The "meta" map, which is a Nextflow/Groovy map with id (the sample name) and a single_end logical entry
    meta = [id: it['sample'], single_end: is_singleEnd]

    // We return a nested map, the first entry is the meta map, the second one is the read(s)
    r2.toString()=='' ? [meta, [r1]] : [meta, [r1, r2]]

}

If you .view() the ch_reads you get:
[[id:sample1, single_end:false], [/Users/atpoint/sample1_R1.fq.gz, /Users/atpoint/sample1_R2.fq.gz]]
[[id:sample2, single_end:true], [/Users/atpoint/sample2_R1.fq.gz]]

A process that can use that could be:

process A {

    tag "$meta.id"

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

    script:

    // The meta map holds the single/paired information that can now easily be parsed in the module in a standardized way
    if(meta.single_end){

        """
        do_single_end_things with ${reads}
        """

    } else {

        """
        do_paired_end_things ${reads}[0] ${reads}[1]
        """

    }

}

A(ch_reads)

The meta map idea (cudos nf-core) is simple, yet versatile and powerful as it allows standardized parsing of such metadata-like information (single_end, libtype, id/sample_name, any other value you can think of) and most importantly, it allows to process single-end paired end data at the same time. All the decision making is based on the meta map and therefore automated if the modules are coded appropriately as above.
Hope that helps.

Welcome to the site @ATpoint


text:
Reference for the statistical methodology used in joint genotype calling
I am searching for a good reference that explains the statistical methodology used in what is variously known as Joint Genotyping, Cohort Genotyping or Multi-sample Genotyping. One with a pseudo code implementation would also be ideal.
I have found several discussing the approach in very general terms and the relative benefits of the technique but nothing with actual equations or algorithms. Failing that, if there is an open source implementation that people would consider to be a reference implementation and that has excellent code documentation that would also be great.

To the best of my knowledge, the first joint calling model was published by Le & Durbin (2011) and refined in Li (2012). They come with detailed derivation. Both samtools and older GATK are using this model. GATK now uses a different model which only works with high-coverage data but is faster to compute. There are other alternatives (e.g. in freebayes).

Here is a recent example of the derivation of the genotype model from Octopus genotype calling paper:

All samples have known ploidy and copy number, so the likelihood function of reads $R$ given genotype $g$ is:
$p(R|g) = \mathop {\prod }\limits_{n = 1}^{|R|} \frac{1}{{|g|}}\mathop {\sum }\limits_{i = 1}^{|g|} p(r_n|g_i)$
where $|g|$ is the ploidy and $|R|$ is the number of reads. The joint genotype posterior for $S$ samples is therefore given by:
$p({\bf{g}}|{{R}},{\cal{M}}_{{g}}) \propto p({\bf{g}}|{\cal{M}}_g)\mathop {\prod }\limits_{s = 1}^S p({{R}}_s|{\bf{g}}_s)$
where the genotype prior model, $\cal{M}$, is either the uniform or HWE-coalescent prior. Unfortunately, the number of genotype combinations g grows exponentially in the number of samples S, so we cannot evaluate the full posterior distribution in general. Therefore, other than for trivial cases, we first approximate the sample marginal genotype posterior distribution under the HWE model (without mutations) $(|,\cal{M})$ and use these marginal probabilities to select K genotype combinations $1,…,$ (K is user-defined) to evaluate under the full joint genotype model. The approximate posterior marginals are computed with expectation maximization.

The source code can be found here.


text:
Maintaining SNP names when calculating allele frequencies from a VCF file on plink 1.9
I downloaded the following vcf file from https://www.internationalgenome.org/.
http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.chr22.phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.vcf.gz
I wanted to first try and calculate allele frequencies, and then from there calculate allele frequencies based on populations and so on. However, I ran into an issue where while I was using plink, I would utilize the --vcf command to try and calculate frequencies. The frequencies were calculated, but the SNP names were missing. I used the following command in Git Bash on Windows:
    plink --vcf ALL.chr22.phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.vcf.gz --freq --out chr22datanew

I then tried to convert the file to plink binary (bed, bim, fam) to try and use that for the basis of the allele frequencies. I used the following commands:
plink --vcf ALL.chr22.phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.vcf.vcf --make-bed --out chr22data

plink --bfile chr22data --freq --out chr22freqs

Once again, the frequencies were calculated with no issue, but the SNP names were missing. I checked the .fam file that was created, and the SNPs were there, but not in the .bim file. This may be a simple, elementary level question, but does anyone know how to fix this? I am using plink 1.9 instead of vcftools because I'm using a computer that runs on Windows OS.

Can you please show the result of head chr22data.bim

Result:
22      .       0       16050075        G       A
22      .       0       16050115        A       G
22      .       0       16050213        T       C
22      .       0       16050319        T       C
22      .       0       16050527        A       C
22      .       0       16050568        A       C
22      .       0       16050607        A       G
22      .       0       16050627        T       G
22      .       0       16050646        T       G
22      .       0       16050654           A

The standard and most common way to name SNPs is by using rsIDs.
The reason your output from PLINK does not include SNP names is because these are not included in the original VCF you are working with. There are at least two possible options in order to fix this.

Annotate first the source VCF with rsIDs from dbSNP. You can do this using bcftools and a recent (or latest) human rsID compilation in VCF format. E.g.

Downloading rsIDs

wget ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh37p13/VCF/All_20180423.vcf.gz
wget ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh37p13/VCF/All_20180423.vcf.gz.tbi

Subsetting chr22 from annotations

This makes the annotation of the chr22 VCF faster
bcftools view --threads 10 --regions 22 All_20180423.vcf.gz -Oz -o 22_20180423.vcf.gz
bcftools index 22_20180423.vcf.gz

VCF Annotation

bcftools annotate --threads 10 -a 22_20180423.vcf.gz -c ID -o ALL.chr22.phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes_annotated.vcf.gz

PLINK frequency estimations

Now you can proceed as you had done already
plink --vcf ALL.chr22.phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes_annotated.vcf.gz --make-bed --out chr22data
plink --bfile chr22data --freq --out chr22freqs

An alternative option is to build your own 'SNP names' based on chromosome and cordinate (e.g. 22_16255382) using plink commands.


text:
umap failed to cluster the cells
I tried umap visualization with scanpy :
sc.pp.normalize_total(adata, target_sum=1e6)
sc.pp.log1p(adata, base=2)
sc.pp.highly_variable_genes(adata, flavor='cell_ranger', n_top_genes=400)

sc.pp.scale(adata, zero_center=True, max_value=None, copy=False, layer=None, obsm=None)
sc.pp.pca(adata, n_comps=50, use_highly_variable=True, svd_solver='arpack')
sc.pp.neighbors(adata, n_neighbors=50)

sc.tl.umap(adata, min_dist=0.5, spread=1.0)
sc.pl.umap(adata, color='fullname', use_raw=False, save='samples_umap.pdf')

version
anndata     0.7.5
scanpy      1.6.1

But the cells can't separate well

I tried another small dataset with scanpy using the same parameters as before:

sc.tl.umap still failed to down dimension the data properly.
Then I tried the original umap package using the same data set:
import umap
import umap.plot
mapper = umap.UMAP().fit(adata.X)
umap.plot.points(mapper)

Now the original umap package can do down dimension very well:

I think there may be something wrong with the umap function in scanpy
Can anyone please let me know the reason?
Thanks a lot.

Without more code and data, this can't be debugged. But first question would be does the data pass basic QC (e.g. look at cellranger web summaries or equivalent). If yes, you could trying different normalisation (code not shown) and/or no scaling. Try sc.tl.umap() with default params. Plot the PCA, try t-SNE. Plot the basic QC metrics following the introductory scanpy tutorial. That is how I would start debugging

I compare the umap in scanpy with the original umap (https://umap-learn.readthedocs.io/en/latest/plotting.html)  using the same dataset, the original umap works well, I think the problem is in scanpy. I edited my question to include this.
Do you have some suggestions?
Thanks

@DanLi It might be this: https://github.com/scverse/scanpy/issues/2045 or related. Try with the latest scanpy version and consider reporting to scanpy github if not already reported

Yes, thanks. @Chris_Rands


text:
VCF: false insertions?
I'm working with gvcf files, containing also non-variant positions. I have noticed some insertions that I'm not able to correctly interpret or handle. An example (selected samples from a huge vcf):
NC_013991.2 1403    .   G   GAA 176,32  .   AC=4;AF=0.014;AN=286;DP=1294;ExcessHet=0.0001;FS=0;InbreedingCoeff=0.5011;MLEAC=3;MLEAF=0.01;MQ=28.09;QD=28.2;SOR=3.912 GT:AD:DP:GQ:PGT:PID:PL:PS   1|1:0,2:2:6:1|1:1403_G_GAA:90,6,0:1403  0/0:2,0:2:6:.:.:0,6,88:.    ./.:0,0:1:.:.:.:0,0,0:. ./.:0,0:0:.:.:.:0,0,0:. ./.:1,0:1:.:.:.:0,0,0:. 0/0:2,0:2:3:.:.:0,3,45:.    ./.:1,0:1:.:.:.:0,0,0:. ./.:1,0:1:.:.:.:0,0,0:. 1/1:0,3:3:9:.:.:123,9,0:.   0/0:1,0:1:3:.:.:0,3,16:.    0/0:1,0:1:3:.:.:0,3,45:.    0/0:2,0:2:6:.:.:0,6,90:.    0/0:9,0:9:21:.:.:0,21,315:.
NC_013991.2 1404    .   A   .   .   .   DP=1286 GT:AD:DP:RGQ    0/0:2:2:6   0/0:2:2:6   0/0:2:2:6   ./.:0:0:0   0/0:1:1:3   0/0:2:2:6   0/0:1:1:3   0/0:1:1:3   0/0:3:3:9   0/0:1:1:3   0/0:1:1:3   0/0:2:2:6   0/0:9:9:21
NC_013991.2 1405    .   A   .   .   .   DP=1282 GT:AD:DP:RGQ    0/0:2:2:6   0/0:2:2:6   0/0:2:2:6   ./.:0:0:0   0/0:1:1:3   0/0:2:2:6   0/0:1:1:3   0/0:1:1:3   0/0:3:3:9   0/0:1:1:3   0/0:1:1:3   0/0:2:2:6   0/0:9:9:21

Here, the reference sequence from position 1403 to 1405 is "GAA". At site 1403, an alternative allele "GAA" is indicated, which I would interpret as an insertion of "AA". Two individuals thus would have genotype "GAAAA" here while the rest has "GAA".
I've seen multiple cases like this, here with longer alternative allele (again, selected samples, not necessarily the same as before):
NC_013991.2 23  .   T   TAATACTTATGTGTTG    66,98   .   AC=2;AF=0.009346;AN=214;DP=314;ExcessHet=0.0105;FS=0;InbreedingCoeff=0.3797;MLEAC=3;MLEAF=0.014;MQ=51.77;QD=33.49;SOR=2.303 GT:AD:DP:GQ:PL  1/1:0,2:2:6:85,6,0  0/0:3,0:3:9:0,9,135 0/0:3,0:3:9:0,9,95  0/0:1,0:1:3:0,3,45  ./.:1,1:2:.:0,0,0   0/0:3,1:4:0:0,0,102
NC_013991.2 24  .   A   .   .   .   DP=322  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:2:2:6   0/0:4:4:12
NC_013991.2 25  .   A   .   .   .   DP=339  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   0/0:4:4:12
NC_013991.2 26  .   T   .   .   .   DP=353  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   ./.:4:5:0
NC_013991.2 27  .   A   .   .   .   DP=361  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   0/0:5:5:15
NC_013991.2 28  .   C   .   .   .   DP=368  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   0/0:5:5:15
NC_013991.2 29  .   T   .   .   .   DP=386  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   ./.:4:5:0
NC_013991.2 30  .   T   .   .   .   DP=388  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   0/0:5:5:12
NC_013991.2 31  .   A   .   .   .   DP=398  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   0/0:5:5:12
NC_013991.2 32  .   T   .   .   .   DP=405  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   0/0:5:5:12
NC_013991.2 33  .   G   .   .   .   DP=419  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   ./.:4:5:0
NC_013991.2 34  .   T   .   .   .   DP=427  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   0/0:5:5:12
NC_013991.2 35  .   G   .   0,09    LowQual BaseQRankSum=-1;DP=433;ExcessHet=3;MLEAC=.;MLEAF=.;MQ=51.42;MQRankSum=1;ReadPosRankSum=0    GT:DP:RGQ   0/0:3:9 0/0:3:9 0/0:3:9 0/0:1:3 0/0:3:9 0/0:5:0
NC_013991.2 36  .   T   .   .   .   DP=441  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   0/0:5:5:12
NC_013991.2 37  .   T   .   .   .   DP=447  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:2:2:6   0/0:3:3:9   0/0:5:5:12
NC_013991.2 38  .   G   .   .   .   DP=449  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:2:2:6   0/0:3:3:9   ./.:4:5:0

The alternative allele "TAATACTTATGTGTTG" is exactly the reference. Of course, the sample in question is homozygous. It seems to me these are "false" insertions: they look like insertions, but aren't.
I can write a simple awk script trimming the the alternative allele, but there are cases where there are real insertions, and I don't know how to distinguish those.
I've tried the solution mentioned in this answer (double pass through bcftools norm) but it doesn't resolve these issues.
I could spend a few days coding a C++ program to do this, but surely someone with better understanding of vcf, C++, etc, has already done this?
Here's a case I really don't know how to interpret:
NC_013991.2 3321    .   A   ATTTCATT    638,04  .   AC=14;AF=0.044;AN=320;DP=1873;ExcessHet=0;FS=0;InbreedingCoeff=0.5521;MLEAC=10;MLEAF=0.031;MQ=29.5;QD=34.42;SOR=1.739   GT:AD:DP:GQ:PGT:PID:PL:PS   1|1:0,2:2:6:1|1:3321_A_ATTTCATT:90,6,0:3321 1|1:0,4:4:12:1|1:3321_A_ATTTCATT:180,12,0:3321  1|1:0,2:2:6:1|1:3321_A_ATTTCATT:90,6,0:3321 0/0:2,0:2:3:.:.:0,3,45:.    0/0:1,0:1:3:.:.:0,3,45:.    0/0:1,0:1:3:.:.:0,3,45:.    0/0:1,0:1:3:.:.:0,3,16:.    0/0:1,0:1:3:.:.:0,3,25:.    0/0:8,0:8:21:.:.:0,21,315:. 0/0:33,0:33:96:.:.:0,96,1440:.  0/0:15,0:15:39:.:.:0,39,585:.
NC_013991.2 3322    .   T   .   .   .   DP=1870 GT:AD:DP:RGQ    0/0:2:2:6   0/0:4:4:12  0/0:2:2:6   0/0:2:2:3   0/0:1:1:3   0/0:1:1:3   0/0:1:1:3   0/0:1:1:3   0/0:8:8:21  0/0:33:33:96    0/0:14:14:36
NC_013991.2 3323    .   T   .   .   .   DP=1859 GT:AD:DP:RGQ    0/0:2:2:6   0/0:4:4:12  0/0:2:2:6   ./.:1:2:0   0/0:1:1:3   0/0:1:1:3   0/0:1:1:3   0/0:1:1:3   0/0:8:8:21  0/0:34:34:99    0/0:14:14:36
NC_013991.2 3324    .   T   .   .   .   DP=1849 GT:AD:DP:RGQ    0/0:2:2:6   0/0:4:4:12  0/0:2:2:6   0/0:1:1:3   0/0:1:1:3   0/0:1:1:3   0/0:1:1:3   0/0:1:1:3   0/0:8:8:21  0/0:34:34:99    0/0:14:14:33
NC_013991.2 3325    .   TGG T   639,64  .   AC=14;AF=0.044;AN=318;DP=1843;ExcessHet=0;FS=0;InbreedingCoeff=0.5416;MLEAC=11;MLEAF=0.035;MQ=29.53;QD=28.53;SOR=1.739  GT:AD:DP:GQ:PGT:PID:PL:PS   1|1:0,2:2:6:1|1:3321_A_ATTTCATT:90,6,0:3321 1|1:0,4:4:12:1|1:3321_A_ATTTCATT:180,12,0:3321  1|1:0,2:2:6:1|1:3321_A_ATTTCATT:90,6,0:3321 0/0:1,0:1:3:.:.:0,3,45:.    0/0:1,0:1:3:.:.:0,3,45:.    0/0:1,0:1:3:.:.:0,3,45:.    0/0:1,0:1:3:.:.:0,3,45:.    0/0:1,0:1:3:.:.:0,3,25:.    0/0:8,0:8:21:.:.:0,21,315:. 0/0:34,0:34:99:.:.:0,99,1485:.  0/0:14,0:14:33:.:.:0,33,495:.
NC_013991.2 3326    .   G   T   1294,06 .   AC=15;AF=0.055;AN=274;BaseQRankSum=0;DP=1833;ExcessHet=0;FS=2.451;InbreedingCoeff=0.5364;MLEAC=17;MLEAF=0.062;MQ=56.66;MQRankSum=0;QD=16.38;ReadPosRankSum=0;SOR=0.828  GT:AD:DP:GQ:PGT:PID:PL:PS   ./.:0,2:2:.:.:.:0,0,0:. ./.:0,4:4:.:.:.:0,0,0:. ./.:0,2:2:.:.:.:0,0,0:. ./.:0,1:1:.:.:.:0,0,0:. ./.:0,1:1:.:.:.:0,0,0:. ./.:0,1:1:.:.:.:0,0,0:. ./.:0,1:1:.:.:.:0,0,0:. 1|1:0,1:1:3:1|1:3309_G_A:45,3,0:3309    1/1:0,4:4:12:.:.:145,12,0:. 0/0:34,0:34:99:.:.:0,99,1485:.  0/0:14,0:14:33:.:.:0,33,495:.
NC_013991.2 3327    .   G   .   0,02    LowQual DP=1832;ExcessHet=3;MLEAC=.;MLEAF=.;MQ=29   GT:DP:RGQ   0/0:2:0 0/0:4:0 0/0:2:0 0/0:1:0 0/0:1:0 0/0:1:0 0/0:1:0 0/0:1:3 0/0:8:18    0/0:34:99   0/0:14:33
NC_013991.2 3328    .   G   A   650,66  .   AC=15;AF=0.048;AN=312;BaseQRankSum=0;DP=1826;ExcessHet=0;FS=3.31;InbreedingCoeff=0.5371;MLEAC=12;MLEAF=0.038;MQ=34.34;MQRankSum=0;QD=27.23;ReadPosRankSum=0;SOR=0.324   GT:AD:DP:GQ:PGT:PID:PL:PS   1|1:0,2:2:6:1|1:3321_A_ATTTCATT:90,6,0:3321 1|1:0,4:4:12:1|1:3321_A_ATTTCATT:180,12,0:3321  1|1:0,2:2:6:1|1:3321_A_ATTTCATT:90,6,0:3321 ./.:0,1:1:.:.:.:0,0,0:. ./.:0,1:1:.:.:.:0,0,0:. ./.:0,1:1:.:.:.:0,0,0:. 0/0:1,0:1:3:.:.:0,3,45:.    0/0:1,0:1:3:.:.:0,3,25:.    0/0:8,0:8:18:.:.:0,18,270:. 0/0:34,0:34:99:.:.:0,99,1485:.  0/0:14,0:14:33:.:.:0,33,495:.
NC_013991.2 3329    .   G   .   18,49   LowQual BaseQRankSum=-1;DP=1818;ExcessHet=3;MLEAC=.;MLEAF=.;MQ=57.63;MQRankSum=0;ReadPosRankSum=1   GT:DP:RGQ   0/0:2:6 0/0:4:12    0/0:2:6 0/0:1:3 0/0:1:3 0/0:1:3 0/0:1:3 0/0:1:3 0/0:7:18    0/0:35:99   0/0:13:33

It seems like only the "CATT" part of the alternative allele is really inserted...
The mapping was done either with bwa aln or bwa mem, and we've observed this behaviour both when using GATK4 and bcftools for genotyping.

"It seems to me these are "false" insertions: they look like insertions, but aren't.": could you explain why they seem false to you? And if we don't need to know about all of the genotypes, please only show us the relevant ones. You say "of course the sample in question is homozygous", but we can't see that since there are so many samples shown. Can you focus only on the relevant ones so we can understand better?

They look false to me because they exactly recapitulate the reference sequence. I removed samples in the examples and could add a few lines illustrating a longer pseudo-insertion.

Sorry if I'm being dense, this sort of thing is always complicated to understand, despite working with this sort of thing every day. Focusing on your first example, the VCF with three lines, we have an insertion of AA, and then two non-variant positions. So some of your samples have GAAAA and some have GAA, the reference, as you say.  What am I missing? What's wrong with that scenario? Isn't this just a case of a duplication indel?

Thanks for sharing your thoughts. I started with the first example as it is short, but most of these "insertions" are rather long such as the second one. I found it suspicious that the large majority of insertions /exactly/ recapitulates the reference sequence. But maybe that's common. Also, we don't seem to get these when mapping RNAseq data with STAR, but with BWA, we have them also in the coding sequences. If the second example really means there are individuals with TAATACTTATGTGTTG and others with TAATACTTATGTGTTGAATACTTATGTGTTG, I invite you to post an answer that I'll accept.

Microsatellites, minisatellites and other type of repetitive variation are a known feature of genomes.
I agree with previous comments in the sense that what you are observing here might as well be real and not just an artefact caused by real repetitive variation in your genome of study, plus read mapping together with variant calling issues at the site.
The latter might be the case, but before getting rid of these alleles I suggest that you inspect the sites and regions around these alleles more carefully: GATK HaplotypeCaller has an option (-- output-bam or similar) that writes a BAM file with the reads as they were used for variant calling after the tool itself performs local re-alignment of reads around the variants. I would use this tool together with the original BAM file in order to inspect the region visually and try to understand it better. For this you can use samtools tview or IGV. The length of your reads becomes somewhat relevant here, but even if you are sequencing around 100 bp, given the small size of these observed alleles you should be able to draw some more solid insights about the variants/errors and how to handle them.

INDELs tend to occur tandemly. Your observation makes them more likely to be real.


text:
Alignment process in Nextflow not picking up the index files from the index process
I'm using bowtie2 to build a small index. After building the index, I want to pass it to the bowtie2 alignment process. The issue is bowtie2-build outputs multiple index files with the .bt2 extension. When I try to provide these index files as an input to the Alignmnet process, it gives the following error:

Error executing process > 'SAM_FOR_STRAND (1)'
Caused by:   Process SAM_FOR_STRAND (1) terminated with an error
exit status (255)
Command executed:
bowtie2 -t -x index.1.bt2 index.2.bt2 index.3.bt2 index.4.bt2
index.rev.1.bt2 index.rev.2.bt2 -1 SRR5655560_1.atria.fastq -2
SRR5655560_2.atria.fastq -S "SRR5655560.sam"
Command exit status:   255
Command output:   (empty)
Command error:   (ERR): "index.1.bt2" does not exist or is not a
Bowtie 2 index   Exiting now ...
executor >  local (4) [71/c095e4] process > INDEX              [100%]
1 of 1 ✔ [1f/35a1e6] process > GTF2BED            [100%] 1 of 1 ✔
[fa/b1e933] process > TRIM (SRR5655560)  [100%] 1 of 1 ✔ [d6/d47ead]
process > SAM_FOR_STRAND (1) [100%] 1 of 1, failed: 1 ✘ Error
executing process > 'SAM_FOR_STRAND (1)'
Caused by:   Process SAM_FOR_STRAND (1) terminated with an error
exit status (255)
Command executed:
bowtie2 -t -x index.1.bt2 index.2.bt2 index.3.bt2 index.4.bt2
index.rev.1.bt2 index.rev.2.bt2 -1 SRR5655560_1.atria.fastq -2
SRR5655560_2.atria.fastq -S "SRR5655560.sam"

The following is simplifies version of the pipeline:
// Declare syntax version
nextflow.enable.dsl=2

/*
 * pipeline input parameters
 */

params.reads = "$projectDir/data/*_{1,2}.fastq"
params.ref_genome = "$projectDir/data/Glycine_max.Glycine_max_v2.1.dna.toplevel.fa"
params.gtf_file = "$projectDir/data/Glycine_max.Glycine_max_v2.1.55.gtf"
params.outdir = "results"

log.info """\

        R N A S E Q - N F   P I P E L I N E
        ===================================
        reads            : ${params.reads}
        reference genome : ${params.ref_genome}
        gtf              : ${params.gtf_file}
        outdir           : ${params.outdir}
        """
        .stripIndent()

/*
 * define the `index` process that creates a binary index
 * given the reference genome file
 */

process INDEX {

  input:
  path ref_genome

  output:
  path "index.*"

  script:
  """
  bowtie2-build $ref_genome index
  """
}

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

process GTF2BED {

    input:
    path gtf

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

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

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

process TRIM {

    tag { pair_id }

    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 1 \\
        -l 12 \\
        -q 30
    """
}

/*
 * define the `SAM_FOR_STRAND` process that creates a SAM file
 * to be used by `STRAND_CHECK` process to determine
 * the strandedness of each fastq file pair.
 */

process SAM_FOR_STRAND {

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

    output:
    path "${pair_id}.sam"

    script:
    """
    bowtie2 -x "${index.baseName}" -1 ${reads[0]} -2 ${reads[1]} -S ${pair_id}.sam
    """
}

workflow {

    index_ch = INDEX(params.ref_genome)

    gtf2bed_ch=GTF2BED(params.gtf_file)

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

    trim_out = TRIM(read_pairs_ch)

    sam_for_strand_ch = SAM_FOR_STRAND(index_ch, TRIM.out.trimmed_fastqs)

}

I've tried passing the index argument directly to the script in the SAM_FOR_STRAND Alignmnet process, it gave the  "index.1.bt2" does not exist error. Then I tried giving the basename of the index files (as shown in the below). It still does not work, even though the index files are being created properly.
 /*
     * define the `SAM_FOR_STRAND` process that creates a SAM file
     * to be used by `STRAND_CHECK` process to determine
     * the strandedness of each fastq file pair.
     */

    process SAM_FOR_STRAND {

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

        output:
        path "${pair_id}.sam"

        script:
        """
        bowtie2 -x "${index.baseName}" -1 ${reads[0]} -2 ${reads[1]} -S ${pair_id}.sam
        """
    }

Best wishes.

When we run Bowtie2, we just need to make sure we supply the index basename using -x. This is the index filename prefix minus the .X.bt2 extension. A safe assumption is that the first file in the collection of output files (sorted in lexicographical order) has the .1.bt2 suffix. We can then just strip this off to get the required prefix. For example:
process SAM_FOR_STRAND {

    input:
    path index_files, stageAs: "bt2_index/*"
    tuple val(pair_id), path(reads)

    output:
    path "${pair_id}.sam"

    script:
    def idx = index_files[0].getBaseName(2)

    """
    bowtie2 \\
        -x "${idx}" \\
        -1 "${reads[0]}" \\
        -2 "${reads[1]}" \\
        -S "${pair_id}.sam"
    """
}


text:
Search ChebiEntity search by common or synonym name?
I can download information from ChEBI if I know the corresponding ChEBI ID:
from libchebipy import ChebiEntity
chebi_entity = ChebiEntity("15903")

And then look at all the synonyms that the compound has:
chebi_entity.get_names()
> [{'_Name__name': 'beta-D-Glucose', '_Name__typ': 'SYNONYM', '_Name__source': 'KEGG COMPOUND', '_Name__adapted': False, '_Name__language': 'en'},
 {'_Name__name': 'beta-D-glucose', '_Name__typ': 'NAME', '_Name__source': 'UniProt', '_Name__adapted': False, '_Name__language': 'en'},
 {'_Name__name': 'beta-D-glucopyranose', '_Name__typ': 'IUPAC NAME', '_Name__source': 'IUPAC', '_Name__adapted': False, '_Name__language': 'en'},
 {'_Name__name': 'BETA-D-GLUCOSE', '_Name__typ': 'SYNONYM', '_Name__source': 'PDBeChem', '_Name__adapted': False, '_Name__language': 'en'},
 {'_Name__name': 'WURCS=2.0/1,1,0/[a2122h-1b_1-5]/1/', '_Name__typ': 'SYNONYM', '_Name__source': 'GlyTouCan', '_Name__adapted': False, '_Name__language': 'en'}]

Is there a way to get the ChEBI ID with one of  the synonyms?
beta-D-glucose -> 15903

Is there any other python tool or web-api that can achieve this?

How very interesting ...
Yes its easy just download the flat file database here:

https://ftp.ebi.ac.uk/pub/databases/chebi/Flat_file_tab_delimited/names.tsv.gz

WARNING: it's a 5.5M file ... don't click the link unless you really want to download a file this size.
You can simply load it straight into a pandas dataframe from read_csv() and then interrogation it via loc. You could place the header NAME as the index. There are other way in, you can shift the dataframe to a dictionary and interrogate that way, or place it straight it SQL via e.g. sqlite3
Here I just loaded it onto Excel and found the entry you were looking for ...

....

Edited: sorry a pandas iloc() solution wasn't appropriate.


text:
How to subset an SRA file for a single chromosome?
I used prefetch to get the Pacbio reads of chicken from the SRA database. I want to align these reads against a reference genome, but not all the reads. I am only interested in a particular region on Chromosome Z. Is there a way that I can subset the downloaded SRA file?
I am trying to use blast, but don't know where to start. How should I prepare the SRA file from prefetch for the blast? I have currently started fasterq-dump for the downloaded SRA file.
Or should I subset after aligning the entire SRA reads data with the reference genome, and in that case what tools do I use?

After dumping out the reads using fasterq-dump, you'll need to align them first and then extract those that map to your region of interest. I think minimap2 is an excellent choice. If you have HiFi reads, you can just use something like:
minimap2 -ax map-hifi ref.fa pacbio-ccs.fq.gz > aln.sam

Then, using samtools, coordinate-sort and index the sorted alignments:
samtools sort -o aln.sorted.bam aln.sam
samtools index aln.sorted.bam

With an indexed BAM, you can then view the alignments that overlap your region of interest:
samtools view aln.sorted.bam chrZ:1-10000

If you instead write these to a file, you can then convert them using samtools  fasta. I think something like the following should suffice:
samtools view -bo chrZ_1_10000.bam aln.sorted.bam chrZ:1-10000
samtools fasta -0 chrZ_1_10000.fa -o /dev/null chrZ_1_10000.bam

You could then upload the FASTA file to the NCBI Web BLAST server.

"After dumping out the reads using prefetch". Do you mean after fasterq-dump ends? And what will be looking at once that ends, a set of fastq files?

@venkateshwar Typo, sorry! Yes, or likely just a single FASTQ file since these are PacBio long reads.


text:
Retrieve rsIDs from chromosome positions using MRutils::get_rsid_from_position
I have 5609 chromosome positions with both ref and alt alleles. I want to retrieve rsIDs for the same. I tried get_rsid_from_position() of the MRutils package in R, but after some time while fetching data from the internet, the cache memory of my system got full, and it froze the system.
Any other way to resolve this issue?
I have 200 GB allotted for the cache memory.

Are you sure you are not incurring into an infinite loop?

An alternative possibility in R is using the library biomaRt from Bioconductor. Here an example starting from a table with 3 columns: Chr, Start, End. Organism is human and positions are assumed to be based on the reference version GRCh38.
library(biomaRt)
snpMart = useEnsembl(biomart = "snps",
         dataset = "hsapiens_snp")

# assuming your table is tab-separated, with column names
myData <- read.table("myTable.txt", sep="\t", header=True)

coords <- apply(myData, 1, paste, collapse = ":")

# submit the query
getBM(attributes = c('refsnp_id', 'chr_name', 'chrom_start', 'chrom_end', 'allele'),
filters = c('chromosomal_region'),
values = coords,
mart = snpMart)

Source

Thank you for suggestion. I have tried this before but getting time-out error even for 500 SNPs. Error in curl::curl_fetch_memory(url, handle = handle) :
Timeout was reached: [www.ensembl.org:443] Operation timed out after 300001 milliseconds with 0 bytes received

@JyotiSharma please show us the actual code you used. An obvious thing to try is 1) just 1 snp, does that work? ii) batch requests: split your list into groups of 20 or 50 or 100 variants and see if that works.


text:
Examples of machine learning approaches to validate results where ground truth is lacking?
I am currently looking into some of the published methods for deconvolution of spatial transcriptomics data where each spot does not have single-cell resolution. These methods all rely on cell type annotated scRNA-seq data to build a reference of expected expression for each cell population. So in essence, because ground truth is more or less impossible to obtain for complete cell typing of a spatially intact tissue section, another type of data (scRNA) is used instead. This lets you infer something about the spatial organization of the tissue despite no spatial information being present in the training data, but in the end no ground truth is available for validating the results.
I was trying to think of an analogous machine learning problem, biological or not, where something similar is done to deal with unattainable ground truch. Are there any obvious/good examples of something like this? Maybe some inspiration can be drawn for how such results are validated.

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.

Hi, discussion and opinion-based questions are too generic to usefully fit the Stack Exchange format. Please adjust your question to relate to a specific issue that you're having. If you want to ask general questions for discussion, it would be better to do that on reddit

If the problem concerns spatial clustering of cells then the ML algorithm is without question DBSCAN.
The problem with deep learning and CNN (convolution) not "deconvolution", i.e. the MBE paper cited, is:

The model is an ad hoc construction, which can be also a strength;
Under-training, is resolved by adding more layers;
Overtraining, can only be assessed post-analysis, in ML its algorithm specific;
CNN approaches require large datasets for training of 2 million+.

The issue with ML is its performance tapers at around the 1 or 2 million data set size mark, whereas deep learning can (stress can) continue to use massive data sets to enhance performance (thats its strength).
Whether, DBSCAN has received publications in bioinformatics - I'm not sure. However, what is present in the public domain isn't always the cutting edge work: a lesson I've very recently learnt. DBSCAN is very fashionable for spatial problems at present and in context personally would rank it above CNN.
CNN vs ML The thing is most of us - unless you're working for Google - don't have access to data sets in 2 million+ range. This can be circumvented by transfer learning but it is not sufficiently developed outside Google and a few major players in computer vision. Has Google performed CNN on one million+ cell images - I doubt it. If the problem is ecological species identification via remote cameras - Googles CNN image libraries would be just the ticket (lots of dogs, cats and bear photos) but for microscopy studies I seriously doubt it.
Conclusion
What I am 100% certain on is DBSCAN is VERY fashionable right now and certain ML algorithms repeatedly outperform deep learning in certain situations. What I am trying to say is ML is not necessarily just a low data size version of deep learning: in some instances it surpasses it and I think DBSCAN has that potential. The broader challenges of deep learning hinge around 'transfer learning'. If biology manages to resolve that then deep learning will start surpassing machine learning, but firstly it's about sharing, secondly in e.g. RNN this is not technically trivial.

What you can do is pose a detailed problem in a later question and I could explain how training would be used ... generally in complex analytical problems it is difficult to assess precisely how a solution would be implemented.


text:
cromwell files downloaded from a URL have too strict permissions
Cromwell can download files that are URLs (i.e. URLs in the inputs.json).  However, the permissions of these files, once they are downloaded, are very strict (rw------- bmac bmac).  I was wondering if there was a way to control the file permissions of said files; I need them accessible to other users.
Cromwell uses a library called akka.http to download url-files.  As stated in the
cromwell configuration docs:

Cromwell uses akka-http to serve requests. For more advanced
configuration settings, refer to the akka-http documentation.

If I run cromwell with the given files
java -jar ~/cromwell/cromwell-84.jar run fq_count.wdl -i fq_count.url.json
I see the permissions for the url-input file:
ls -l cromwell-executions/fq_count/2a65649a-999e-44ef-bec8-89953e0e3f1f/call-count_reads/inputs/429920399/12673489687364390176
-rw-------  2 bmac  bmac

Question:  Is there a way to control the file permissions, i.e. through akka.http or other?
WDL
version 1.0
workflow fq_count {
    input {
        File fastq_file
    }
    call count_reads {
      input: fastq = fastq_file
    }
}

task count_reads {
  input {
    File fastq
    String memory = "5G"
    String time = "00:10:00"
  }

  command <<<
    set -euo pipefail
    fastq='~{fastq}'
    extension="${fastq##*.}"

    if [ "${extension}" = 'gz' ]; then
  cat_cmd='zcat'
else
  cat_cmd='cat'
fi
# using awk instead of bc as bc isn't in debian slim images
"${cat_cmd}" "~{fastq}" | wc -l | awk '{ print int($1/4) }' > num_reads.txt
  >>>

  output {
    File num_reads = "num_reads.txt"
  }

  runtime {
    cpu: 1
    memory: memory
    time: time
  }
}

Input json
{
  "fq_count.fastq_file": https://portal.nersc.gov/cfs/m342/jaws/test_data/sample.fastq"
}

Is adding a chmod command to the end an option? Ugly, but it might work.


text:
How to get separate histograms plots on the basis of the column value?How to detect which plot has most deviation?
I have a data frame (df) which has correlations calculated for different genes with respect to different ID combinations. I want to get separate histogram plot based on  the gene name (separate plot for separate gene). Also after plotting I want to know which genes have more deviations in their correlation values (Like more deviated from the median in each plot).
Example:
df
Gene    ID_1    ID_2    Correlation
TNF 1175077258  1175075956  0.626115
TNF 1175077261  1175075956  0.542499
TNF 1175077261  1175077258  0.50113
TNF 1175077263  1175075956  0.431587
TNF 1175077263  1175077258  0.734858
TNF 1175077263  1175077261  0.622206
CCL2    952399885   952399878   0.632703
CCL2    952399894   952399878   0.622244
CCL2    952399894   952399885   0.697313
CCL2    952399923   952399878   0.687826
CCL2    952399923   952399885   0.612089
CCL2    952399923   952399894   0.603834
IL10    1068768000  1068767927  0.118955
IL10    1068768147  1068767927  0.32038
IL10    1068768147  1068768000  0.430287
IL10    1068768409  1068767927  0.335264
IL10    1068768409  1068768000  0.406426
IL10    1068768409  1068768147  0.546452

I tried the below code which only works when combining all genes together.
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams["figure.figsize"] = (8,6)

g1 = sns.histplot(df.Correlation)
g1.set(
    xlim=(-1,1),
    title="correlation"
)
g1.axvline(0.0, c='black', linestyle="--")
g1.axvline(df.Correlation.median(), c='red', linestyle="--")

The result/plot I am expecting is this (I have created this example plot in excel)

You could use the groupby() pandas function to group the dataframe by gene name. And then just loop through each group to plot the histogram of correlation values. I think you can use the mean of the ratio (deviation / median) to compare deviation for different genes:
import matplotlib.pyplot as plt

for name, group in df.groupby('Gene'):

    median = group['Correlation'].median()
    deviation = abs(group['Correlation'] - median)
    ratio = deviation / median
    ratio_mean = ratio.mean()

    id_1 = group['ID_1'].astype(str)
    id_2 = group['ID_2'].astype(str)

    group['Combination'] = id_1 + '/' + id_2

    plt.figure(figsize=(8, 6))
    plt.bar(group['Combination'], group['Correlation'])

    plt.xlabel('Combination', fontsize=12, fontweight='bold')
    plt.ylabel('Correlation', fontsize=12, fontweight='bold')

    plt.xticks(rotation=45, ha="right")
    plt.title(f'{name} (deviation ratio: {ratio_mean.round(3)})')
    plt.tight_layout()

    plt.gca().spines['top'].set_visible(False)
    plt.gca().spines['right'].set_visible(False)

    plt.show()

Is there a way to also get correlation matrix based on the above df data frame for each gene separately? Thank you

@Riya I think what you want is the corr() function for this: corr_matrix = group[['Combination', 'Correlation']].corr()

Thank you Steve


text:
dealing with a list of VEP files
I'm new to dealing with big and compressed data so, I have some questions regarding that
I have 1000 files, each file has the extension (vep.txt.gz) I want to read these files, then filter them. The size of these files is very large. I tried to read these files in Python and R but I am facing a memory error issue.
I need to read these files and extract some rows from each file for example I have two columns. In the first column, I want the rows with low values in the impact column and the rows with likely pathogenic values from the clinSIG column. Also, I have to add a new column to put the files' names as IDs for each sample.
So, how can I read a list of big (vep.txt.gz) files in R or Python, then filter the content of the files based on two columns (impact = low) and (clin_SIG = Likely pathogenic) ?
Should I decompress the files or just read them?
The following is an example of the files but each file is compressed and all the files have the same header (you can find the original file, which is tab separated, [here][1]):
## ENSEMBL VARIANT EFFECT PREDICTOR v108.0
## Output produced at 2017-03-21 14:51:27
## Connected to homo_sapiens_core_108_38 on ensembldb.ensembl.org
## Using cache in /homes/user/.vep/homo_sapiens/108_GRCh38
## Using API version 108, DB version 108
## polyphen version 2.2.2
## sift version sift5.2.2
## COSMIC version 78
## ESP version 20141103
## gencode version GENCODE 25
## genebuild version 2014-07
## HGMD-PUBLIC version 20162
## regbuild version 16
## assembly version GRCh38.p7
## ClinVar version 201610
## dbSNP version 147
## Column descriptions:
## Uploaded_variation : Identifier of uploaded variant
## Location : Location of variant in standard coordinate format (chr:start or chr:start-end)
## Allele : The variant allele used to calculate the consequence
## Gene : Stable ID of affected gene
## Feature : Stable ID of feature
## Feature_type : Type of feature - Transcript, RegulatoryFeature or MotifFeature
## Consequence : Consequence type
## cDNA_position : Relative position of base pair in cDNA sequence
## CDS_position : Relative position of base pair in coding sequence
## Protein_position : Relative position of amino acid in protein
## Amino_acids : Reference and variant amino acids
## Codons : Reference and variant codon sequence
## Existing_variation : Identifier(s) of co-located known variants
## Extra column keys:
## IMPACT : Subjective impact classification of consequence type
## DISTANCE : Shortest distance from variant to transcript
## STRAND : Strand of the feature (1/-1)
## FLAGS : Transcript quality flags

#Uploaded_variation  Location   Allele  Gene             Feature          Feature_type  Consequence                            cDNA_position  CDS_position  Protein_position  Amino_acids  Codons   Existing_variation  IMPACT    DISTANCE  STRAND  CLINSIG
11_224088_C/A        11:224088  A       ENSG00000142082  ENST00000525319  Transcript    missense_variant                       742            716           239               S/I          aGc/aTc  -                   LOW  -         -1      likely pathogenic
11_224088_C/A        11:224088  A       ENSG00000142082  ENST00000534381  Transcript    downstream_gene_variant                -              -             -                 -            -        -                   MODIFIER  1674      -1      -
11_224088_C/A        11:224088  A       ENSG00000142082  ENST00000529055  Transcript    downstream_gene_variant                -              -             -                 -            -        -                   HIGH  134       -1      -
11_224585_G/A        11:224585  A       ENSG00000142082  ENST00000529937  Transcript    intron_variant,NMD_transcript_variant  -              -             -                 -            -        -                   MODIFIER  -         -1      -

From the previous file, the expected output will be one row which is the first row because the impact is low & CLINSIG is likely pathogenic so i want to find all rows with this condition for 1000 files
ID  #Uploaded_variation  Location   Allele  Gene    Feature Feature_type  Consequence   cDNA_position   CDS_position    Protein_position  Amino_acids  Codons   Existing_variation  IMPACT    DISTANCE  STRAND  CLINSIG
M1457.vep   11_224088_C/A   11:224088   A   ENSG00000142082 ENST00000525319 Transcript  missense_variant    742 716 239 S/I aGc/aTc -   LOW  -  -1  likely pathogenic

Thanks in advance

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.

What is in vep.txt.gz exactly? How did you try to read this in Python? What, exactly, do you want to filter and how? What are your filtering criteria? If you show us an example input (you can decompress the file and post a few lines here) and the output you want from it, then we can help you but without any information about the file or its format, we cannot really do much.

I've edited the question hope it's now more clear and also hope you can help me @terdon

Thanks but I'm afraid we can't do anything with an image. We need to have an example that we can use to test our solutions. Please add the file as text (paste it into your question, and use the {} button to format it as code) and also show us the output you expect from it so we can check if we get it right.

now you will find the file @terdon

Should we assume the file is tab-separated?

yeah, the format is tab-separated @terdon

Hang on, what output do you want? You say "the expected output will be [. . . ] the first row" but earlier you said "I have to add a new column to put the files' names as IDs for each sample". So what do you actually want as output? Please edit your question and add your expected output.

Let us continue this discussion in chat.

For these requirements, it might be possible to simply use filter_vep included in the VEP package for this. Note that multiple --filter flags are ANDed together (i.e. all filters must pass for a line to be printed). For example:
for i in *.vep.txt.gz; do
    filter_vep \
        -i "${i}" \
        -o "${i%.vep.txt.gz}.filtered.vep.txt" \
        --filter "IMPACT is LOW" \
        --filter "CLIN_SIG matches likely_pathogenic"
done

If you then wanted to merge the output and include the filename, you could just use awk. With the following in a file called merge_filtered_vep_txt.awk:
BEGIN {
    FS=OFS="\t"
    f=1
}

/^##/ {
    next
}

sub(/^#/, "") && f {
    print "ID", $0
    f=0
    next
}

{
    print FILENAME, $0
}

Run using:
awk -f merge_filtered_vep_txt.awk /path/to/*.filtered.vep.txt

Python can also open gzip-compressed files for reading (and writing), using the gzip module in the standard library. These can be parsed line-by-line to avoid memory errors. Here's one way using the csv module (also in the standard library):
import argparse
import contextlib
import csv
import gzip
import pathlib
import sys

def get_argument_parser():

    parser = argparse.ArgumentParser(add_help=False)

    group = parser.add_argument_group('input options')
    group.add_argument(
        "-i",
        "--input",
        dest='glob_pattern',
        type=str,
        default="**/*.vep.txt.gz",
        metavar="STR",
        help="The relative glob pattern (default: '%(default)s')",
    )
    group.add_argument(
        "-d",
        "--dir",
        dest='glob_path',
        type=pathlib.Path,
        metavar="DIR",
        default='.',
        help="Glob files in DIR (default: '%(default)s')",
    )

    group = parser.add_argument_group('output options')
    group.add_argument(
        "-o",
        "--out",
        type=pathlib.Path,
        metavar="FILE",
        default='-',
        help="Write the filtered output to FILE (default: stdout)",
    )
    group.add_argument(
        "-f",
        "--force",
        action='store_true',
        help="Overwrite the output file if it exists",
    )

    group = parser.add_argument_group('additional options')
    group.add_argument(
        "-h",
        "--help",
        action="help",
        help="Show this help message and exit",
    )

    return parser

def main():

    parser = get_argument_parser()
    args = parser.parse_args()

    with contextlib.ExitStack() as stack:
        if args.out.name != '-':
            out = stack.enter_context(open(args.out, 'wt' if args.force else 'xt'))
        else:
            out = sys.stdout

        # Get a list of files from the directory specified by args.glob_path

        file_list = sorted(args.glob_path.glob(args.glob_pattern))

        # Build up a unique list of header fields from each input file

        header_cols = ['ID']
        for gz_file in file_list:
            with gzip.open(gz_file, 'rt') as infile:
                for line in infile:
                    if not line.startswith('##'):
                        break

                assert line.startswith('#')

                for field in line[1:].rstrip().split('\t'):
                    if field not in header_cols:
                        header_cols.append(field)

        # Create a csv.DictWriter() and write the header fields

        writer = csv.DictWriter(out, fieldnames=header_cols, delimiter='\t')
        writer.writeheader()

        # Now re-open each file in the arguments list for reading

        for gz_file in file_list:
            with gzip.open(gz_file, 'rt') as infile:

                # Skip the info header lines
                line = next(infile)
                while line.startswith('##'):
                    line = next(infile)

                # Get the column headers
                fieldnames = line.rstrip().removeprefix('#').split('\t')

                # Parse each input file and store the desired rows in memory
                reader = csv.DictReader(infile, fieldnames=fieldnames, delimiter='\t')
                rows = []
                for row in reader:

                    impact = row.get('IMPACT')
                    clinsig = row.get('CLINSIG')

                    if impact == 'LOW' and clinsig == 'likely pathogenic':
                        rows.append({'ID': gz_file, **row})

                # Write the list of rows to the output file
                writer.writerows(rows)

if __name__ == '__main__':
    main()

With the above in a file called filter_vep_txt.py:
usage: filter_vep_txt.py [-i STR] [-d DIR] [-o FILE] [-f] [-h]

input options:
  -i STR, --input STR  The relative glob pattern (default: '**/*.vep.txt.gz')
  -d DIR, --dir DIR    Glob files in DIR (default: '.')

output options:
  -o FILE, --out FILE  Write the filtered output to FILE (default: stdout)
  -f, --force          Overwrite the output file if it exists

additional options:
  -h, --help           Show this help message and exit

thank you in python code, what is the line code for determining the list of files of vep.txt.gz? @Steve

vep.txt.gz? ----> I suspect is just one file named vep.txt see: https://stackoverflow.com/questions/31028815/how-to-unzip-gz-file-using-python ; https://unix.stackexchange.com/questions/93139/can-i-zip-an-entire-folder-using-gzip

@Maj Maybe I haven't understood correctly, but I think that's really a question for you. I.e. what files would you like to process? How would you like to come up with the list? Would you prefer some sort of glob pattern or would you prefer to supply a file of filenames, one per line? Perhaps you could edit your question to include some details?

I've edited the question hope it's now more clear and also hope you can help me @Steve

@Maj Please see the updated code above. It's not clear if each file contains the same header fields, so the above reads the list of input files twice to ensure we get a unique list of header columns. Run using python filter_vep_txt.py -o out.txt *.vep.txt.gz or similar.

thanks a lot @Steve, sorry but I can not know how your code finds the files. is it for one files or for all?

What is the command to locate the files? @Steve

@Maj You'll need to supply one or more files as positional arguments on the command line. If all of your files are in the current directory, you should be able to just use the above glob pattern. If these are under multiple directories, you could use: python filter_vep_txt.py -o out.txt /path/to/files/*.vep.txt.gz /path/to/more/files/*.vep.txt.gz for example.

ok thank you @Steve now it is clear, I saved the first 4 blocks of python code in filter_vep_txt.py but the last block code you wrote how can I run it? sorry I'm new in this field

@Maj the last block is just the usage. You can print it using: python filter_vep_txt.py -h

@Maj In my code above, you may need to change 'Likely pathogenic' to 'likely pathogenic'. Is your CLINSIG column actually a comma-separated list? I'm just curious, but what options did you use to run VEP? Maybe you could include this in your question?

hi @Steve, when I run the code by python filter_vep_txt.py -o out.txt .vep.txt.gz, It shows me this error>> OSError: [Errno 22] Invalid argument: '.vep.txt.gz'.     NOTE (all the files are in the same directory)

@Maj You just have a typo in your command, use: python filter_vep_txt.py -o out.txt *.vep.txt.gz. Note the star/asterisk in the glob pattern.

@Steve yeah I tried with one file and with all files but still I get the same error ((OSError: [Errno 22] Invalid argument: '*.vep.txt.gz'))

Traceback (most recent call last):
File "C:\Users \Desktop\files\filter_vep_txt.py", line 123, in
main()
File "C:\Users\Desktop\files\filter_vep_txt.py", line 75, in main
with gzip.open(vep_file, 'rt') as infile:
File "C:\Program Files\WindowsApps\PythonSoftwareFoundation.Python.3.10_3.10.2544.0_x64__qbz5n2kfra8p0\lib\gzip.py", line 58, in open     binary_file = GzipFile(filename, gz_mode, compresslevel)
File "C:\Program Files\WindowsApps\PythonSoftwareFoundation.Python.3.10_3.10.2544.0_x64__qbz5n2kfra8p0\lib\gzip.py", line 174, in init
@Steve

fileobj = self.myfileobj = builtins.open(filename, mode or 'rb')
OSError: [Errno 22] Invalid argument: '*.vep.txt.gz'
@Steve sorry for many questions

@Maj I think the problem is just that your shell can't expand the glob pattern. I think the simpler solution is to just have the script find the files. Please see the changes above. Use like: python filter_vep_txt.py -d /path/to/files -o out.txt, and it will try to glob the files recursively.

@Steve the code is working now, but when I open the out.txt file, I found only the word (ID) and I did not find any other information, I ran this command  (python filter_vep_txt.py -d *.vep.txt.gz -o out.txt) in the same files path

@Maj Use -i instead of -d: python filter_vep_txt.py -i '*.vep.txt.gz' -o out.txt. Be sure to use single quotes to prevent it from being expanded by your shell. Failing that, add a print(line) before the assert statement and post the output. I'm not expecting VEP output files to contain blank lines in the header, but maybe your output is slightly different.

@Steve thanks a lot for your reply... I get this error and when I open out.txt file I found only the header of files (columns names) ..... Traceback (most recent call last):
File "C:\Users\t\Desktop\files\filter_vep_txt.py", line 122, in
main()
File "C:\Users\t\Desktop\files\filter_vep_txt.py", line 116, in main
rows.append({'ID': vep_file, **row})
NameError: name 'vep_file' is not defined. Did you mean: 'gz_file'?

@Maj Sorry, typo. Yes I did indeed mean 'gz_file'.

You don't need R or Python for this, it is a very straightforward text parsing problem and you can easily do it in an awk one-liner:
$ awk 'BEGIN{FS=OFS="\t"}
      /^[^#]/ && $14=="LOW" && $17="likely pathogenic" ' file
11_224088_C/A   11:224088   A   ENSG00000142082 ENST00000525319 Transcript  missense_variant    742 716 239 S/I aGc/aTc -   LOW -   -1  likely pathogenic

Explanation

BEGIN{FS=OFS="\t"}: the BEGIN{} block is run before processing the file. Here we are setting the field separator (FS) and output field separator (OFS) to a tab so that awk will split each input line on tab characters.
/^[^#]/ && $14=="LOW" && $17="likely pathogenic": this will evaluate to true if the current line does not start with a #, and the 14th field is LOW and the 17th field is likely pathogenic. In awk, the default action when something evaluates to true is to print the current line, so this results in the expected output.

Now, if you also want to add the file name as a field, you can do this to have it as the last field:
$ awk 'BEGIN{FS=OFS="\t"} /^[^#]/ && $14=="LOW" && $17="likely pathogenic"{$(NF+1)=FILENAME; print} ' file
11_224088_C/A   11:224088   A   ENSG00000142082 ENST00000525319 Transcript  missense_variant    742 716 239 S/I aGc/aTc -   LOW -   -1  likely pathogenic       yourfilename

Next, if your files are compressed and you want to do this to all of them, you could simply do:
for file in *vep.txt.gz; do
  zcat "$file" | awk -v fname="${file%.vep.txt.gz}" '
                       BEGIN{FS=OFS="\t"}
                       /^[^#]/ && $14=="LOW" && $17="likely pathogenic"{
                          $(NF+1)=FILENAME;
                          print
                       }'
done > allfiles.filtered

I would just use the included filter_vep script and merge the output files using AWK if needed. I'm not sure if VEP always produces it's output fields in exactly the same order, so using an array to refer to the IMPACT and CLINSIG fields might be preferable.

@Steve sure, but you already have given a good answer using the included tool, and I personally find it easier and more intuitive to use general scripting tools instead of taking the time to learn tool-specific helper scripts each time. Fair point about the array, but I admit I would be shocked if VEP is so badly written that it changes the order of its fields!

Assuming it's been run the same way for each file you're probably fine. But if you were to re-run it with a different set of options, the order of the non-default fields (incl. IMPACT and CLIN_SIG) could change. For me, using a recent VEP, the non-default fields actually end up in an "Extra" column consisting of semicolon-separated key=value pairs.

@terdon do you know, how I can run these commands in Windows 11?   for file in *vep.txt.gz; do...etc

No, @maj, sorry I don't use Windows. You should be able to do it using windows subshell for windows.

@Maj I think terdon means the Windows Subsystem for Linux, WSL2

I did indeed, thanks @Steve

thank you @terdon  The whole code worked correctly except for the command to add a column for the filename as an identifier, do you know what is the problem?

@Maj sorry but I am not sure what you mean. If you are referring to the FILENAME part, it is possibel that you are using an awk that doesn't support it. GNU awk (gawk), which is the default on Linux systems, supports it but I'm not sure about others.


text:
How to get a file with the introns of a gene?
I was able to get the exons of a particular gene as a fasta file from the NCBI genome viewer. As for introns I only found a bed file that showed where the intron starts and ends, along with its direction. What tool do I use to get the exact sequence of these 11 introns using this bed file?
And is there a more direct way?
Thank you.

First I edited the Chromosome column of the bed file from ncbi, because it contained the accession number instead of the chromosome. Then I used the getfasta tool of the bedtools utility.
bedtools getfasta -fi input_genome_file.fa -fo output.fa -bed bed_file
The output contained the introns in the format I wanted.
If anybody had a better way, feel free to post.

Hi @venkateshwar thank you for sharing this information. Could you just confirm that this is the answer to the question you posted (it looks like it)?

@M__ Ok I will.