text:
How deep variant and others preprocessed insertion and deletion for extracting candidate variant?
I am trying to build a model for variant calling,but I am a bit confused how to pre-process the reads that has insertion and deletion and align them with reference?
Should I add something to the location of insertion or deletion to reference or should I just delete from read?

Can you provide some reference for why you think preprocessing is necessary?

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.

Are you building a model specifically for deep variant, or trying to create your own model via a different algorithm?


text:
Installation of Biovia DS visualiser on ubuntu
I am trying to install Biovia DS visualiser2021 on ubuntu 22.04. I used the procedure from the following website
https://askubuntu.com/questions/1175772/cant-install-discovery-studio-visualizer
When I run the installer ./install_DSClient.sh I got the following error
Calling LicensePack Installer

Verifying archive integrity... All good.

Uncompressing License Pack and extracting files, please wait...
./install_lp.sh: 73: Syntax error: redirection unexpected
Copying 1 files to ./share/license/data/tempdata.dat Please wait...
Copying 1 files to /home/user/BIOVIA/DiscoveryStudio2021/Logs/dsclientinstall.log Please wait...
Completed installing  DiscoveryStudio2021 at /home/user/BIOVIA/DiscoveryStudio2021

Check log file [ /home/user/BIOVIA/DiscoveryStudio2021/Logs/dsclientinstall.log ] for more information. Return code: 0

How to solve this error?

source: Stack Exchange

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

sort: votes


text:
Multiple sequence alignment with gaps
I'm working on a problem where the elements of multiple arrays need to be aligned with each other. While doing some research on algorithms that could possibly be applicable to this task, I stumbled on sequence alignment. I believe this is very similar to what I'm trying to achieve.
In short, I need to have equal elements aligned with each other across multiple arrays. However, the difference with my problem compared to sequence alignment seems to be that the latter doesn't tolerate gaps. On the other hand, gaps would be perfectly acceptable for me. So, for example, given input:
[A, ,B,A, ]
[ ,C, ,D,A]
[B, ,C,B,A]
A desired output would be:
[A,B, , , ,A]
[ , ,C, ,D,A]
[ ,B,C,B, ,A]
But there cannot be any errors. That is, every column must contain elements of the same value.
Would it be possible to achieve the desired result with sequence alignment algorithms?

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


text:
What statistical analysis should I perform on DNA variants from unequal number of case (34) and control (11) samples?
I have performed NGS on 34 patient samples and 11 controls. Now I have variants from both groups but there is obvious difference in number of variants.
How can I compare variants of both groups statistically?
Can I do association analysis?

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.

This is simply a GWAS calculation. You can do this via phylogenetics, it's very complicated, but it will work. HOWEVER, I am assuming this is a whole genome NGS (variants so it must be), rather than RNAseq. RNAseq is a very different calculation. It would be useful to specify the disease because there will be prior variant association studies.

Please explain more about your experiment and desired analysis. What platform was sequencing done on? What coverage? How do you know you have variants in both groups? What additional processing has been done on the FASTQ files to demonstrate that there are variants?

Mitochondrial DNA sequencing on Illumina platform and analysis was performed using bwa-mem, bcftools and Mitomap for variant annotation. I want to compare DNA variants in both groups (34 patients, 11 controls) and any association of these variants with disease. Can population specific Allele frequency from gnomAD help in analysis of rare variants?


text:
Identification of unique chains in proteins by sequence similarity
I'm trying to find unique chains or proteins within their PDB files. Many proteins have multiple chains, but very often they are identical (say, the PDB file consists of 4 homodimers, for example, so there is one unique chain represented 8 times total). I am processing certain structural features of the chains, so processing two identical chains is redundant for my pruposes.
I am using Biopython and the way I'm currently identifying unique chains is looking through each PDB file, checking the sequences of the chains and comparing them to one another using alignments (specifically pairwise2.align.globalxx(seq1, seq2). I consider a chain non-unique if it shares >90% sequence similarity to another chain. (Most chains have 100+ residues.)
My question is, is this threshold appropriate? Should I go lower, maybe 80%, or the opposite, higher? For now, it seems to work fine, but I'll be processing a lot of PDB files and want to be on the safer side. If I accidentally process identical chains, that's not much of an issue, but it adds unnecessary time to the processing. Is there maybe a source/article to help with what sequence similarity I should be using? I've found multiple articles, but they usually talk about homology (i.e., common ancestry of proteins), but I'm only doing this within the protein itself and its subunits (chains), so this is of no use for me.

I commented yesterday w/ a link to SIFTS —sorry about the terseness: I meant to say why not use SIFTS table of PDB chain to Uniprot if you are just after finding which chains are the same. Alternatively, the PDB API data will tell you in a more convoluted way what is present. Neither will tell you about gaps in sequence and the latter will be weird for C-Alpha traces (as the precise AA is unknown)

Came across this looking for a solution myself last night but now I have created my own I will share it.
Rather then using Biopython (which is a Python package I use a lot), it is more computationally efficient to parse the .fasta files from the protein data bank (RCSB in this case) which as far as I can see, contain this information. For example, this is 2P16's .fasta file (with 2 unique chains):
>2P16_1|Chain A|Coagulation factor X (EC 3.4.21.6) (Stuart factor) (Stuart-Prower factor)|Homo sapiens (9606)
IVGGQECKDGECPWQALLINEENEGFCGGTILSEFYILTAAHCLYQAKRFKVRVGDRNTEQEEGGEAVHEVEVVIKHNRFTKETYDFDIAVLRLKTPITFRMNVAPACLPERDWAESTLMTQKTGIVSGFGRTHEKGRQSTRLKMLEVPYVDRNSCKLSSSFIITQNMFCAGYDTKQEDACQGDSGGPHVTRFKDTYFVTGIVSWGEGCARKGKYGIYTKVTAFLKWIDRSMKT
>2P16_2|Chain B[auth L]|Coagulation factor X (EC 3.4.21.6) (Stuart factor) (Stuart-Prower factor)|Homo sapiens (9606)
KLCSLDNGDCDQFCHEEQNSVVCSCARGYTLADNGKACIPTGPYPCGKQTLE

and 6WON's (with 8 of the same chains):
>6WON_1|Chains A, B, C, D, E, F, G, H|YohF|Salmonella typhimurium (strain LT2 / SGSC1412 / ATCC 700720) (99287)
QGMTKVAIVTASDSGIGKACALLLAQNGFDIGITWHSDERGAQETAKKAAQFGVRAETIHLDLSQLPEGAQAIEHLIQRLGRVDVLVNNAGAMTKSAFIDMPFTQWRQIFTVDVDGAFLCAQIAARHMIKQGEGGRIINITSVHEHTPLPQASAYTAAKHALGGLTKSMALELIEHHILVNAVAPGAIATPMNDMDDSDIEPGSEPSIPIARPGSTHEIASLVAWLCSEGASYTTGQSLIVDGGFMLANPQFNAK

So far this always seems to be the case and parsing the .fasta is a lot easier then using BLAST or sequence similarity such as you are or were.
For the implementation, whereas the .pdb, .cif files etc use the same download URL (https://files.rcsb.org/download/), the .fasta files use this link: https://www.rcsb.org/fasta/entry/. Parsing the downloaded .fasta files and only reading the lines that begin with >, stripping out unnecessary chains (in my specific case I just want the first unique chains in the PDB so A, B being the same and C, D being the same would return me a list with [A, C], you can modify this script to your own ends) and running the Python code gives this output:
First unique chains in 1d2z (2 unique chains): ['B', 'A']
First unique chains in 2p16 (2 unique chains): ['A', 'B']
First unique chains in 6won (1 unique chain):  ['A']

There was an issue with some entries having Chain or Chains so both had to be considered. Additionally, there are commas that can be removed from all but the last chain ID but for my purposes I do not need to strip it, you can put something in _fasta_helper() to this end before return. This is the Python code and hopefully it helps someone out in the future:
import re
from typing import Union
import urllib.request
from urllib.error import (HTTPError,
                          URLError,
                          ContentTooShortError as ContentError)

def rcsb_downloader(pdb: str, ext: str = None, down_path: str = None) -> str:
    """Downloads the PDB from the protein databank RCSB.

    Parameters
    ----------
    pdb : str
        Protein databank accession 4 letter code i.e. '6won'
    ext : str, default = None
        Extension of the file to download i.e. '.fasta', '.cif' etc. Note: if
        None argument then ext = '.pdb'

    Returns
    -------
    str | bool
        Returns str path if download successful, False if failure
    """

    if ext is None:
        ext = '.pdb'
    if down_path is None:
        download_path = ('/your/path/to/dir/' + pdb + ext)
    else:
        download_path = down_path + pdb + ext

    rcsb_path = 'https://files.rcsb.org/download/' + pdb + ext

    if ext == '.fasta': # fasta has a different web search then others
        rcsb_path = 'https://www.rcsb.org/fasta/entry/' + pdb

    try:
        urllib.request.urlretrieve(rcsb_path, download_path)
        return download_path
    # Try again in case there was a temporary issue with the first attempt
    except (ContentError, HTTPError, URLError):
        try:
            urllib.request.urlretrieve(rcsb_path, download_path)
            return download_path
        except (ContentError, HTTPError, URLError):
            print(f"PDB failed to download {pdb}")
            return False

    def _fasta_helper(self, line: str, single_unique_chain: bool = True) -> Union[str, list]:
        """Separate the first chain ID from the .fasta comment line."""

        # [1] some of the chains had [auth C] etc after so this is removed
        # [2] was originally a list comprehension that filtered based on the word
        # chain, however some of the fasta entries had "Chain" or "Chains" in their
        # description which broke this. The first arrow is removed, split at |, the
        # second entry (1st in Python) entry contains the chain information and the
        # first entry removed (which so far has always been chain, anticipate their
        # being edge cases with older PDB entries). This appears to hold up.

        clean_chain = [
            x.replace(",", "")
            for x in list(
                re.sub(r"\[.*?\]", "", line).lstrip(">").split("|")[1].split(" ")[1:]  # [1]  # [2]
            )
        ]

        if not single_unique_chain:
            return clean_chain
        return clean_chain[0]

def fasta_parse(path_: str) -> list:
    """Parse .fasta files and analyse lines beginning with '>'."""

    with open(path_, 'r', encoding='utf-8') as fasta:
        total_unique_chains = 0  # Use as a flag to test if more than one unique chain
        unique_chains       = []
        for line in fasta:
            if line[0] == '>' and total_unique_chains == 0:
                total_unique_chains += 1
                unique_chains.append(_fasta_helper(line))
                continue
            elif line[0] == '>' and total_unique_chains != 0:
                unique_chains.append(_fasta_helper(line))
    return unique_chains

pdb_diff_chain  = '2p16' # Chain ...
pdb_diff_chains = '1d2z' # Chains ...
pdb_same_chains = '6won' # Chains ...

download_p_diff  = rcsb_downloader(pdb_diff_chain, '.fasta')
download_p_diffs = rcsb_downloader(pdb_diff_chains, '.fasta')
download_p_same  = rcsb_downloader(pdb_same_chains, '.fasta')

unique_s  = fasta_parse(download_p_same)
unique_d  = fasta_parse(download_p_diffs)
unique_ds = fasta_parse(download_p_diff)

print(f'First unique chains in {pdb_diff_chains} (2 unique chains): {unique_ds}\n'
      f'First unique chains in {pdb_diff_chain} (2 unique chains): {unique_d}\n'
      f'First unique chains in {pdb_same_chains} (1 unique chain):  {unique_s}')

This is a great way, but unfortunately I have the files as PDBs, so I'd have to re-download everything as FASTA files instead. However, I really like your solution, thank you for it! It's a great way to do it and definitely better than processing sequences from PDBs (due to chain breakings and its processing by BioPython.Polypeptide module), so I'll keep this in mind when I need to do it again.


text:
How can I download from NCBI all the ITS genes and the related taxonomy?
I would like to download all the ITS1 and ITS2 genes from NCBI in a fasta file. And, I'd like to download even the related taxonomy of each sequence.
Thanks,
Marco

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

I recommend to use Biopython:
from Bio import Entrez

Entrez.email = "your@somemail.com"
Entrez.api_key = "YOUR_API_KEY"

search_result = Entrez.read(Entrez.esearch(db="nucleotide", term="ITS1[Gene Name]"))
with Entrez.efetch(
    db="nucleotide", id=search_result["IdList"], rettype="fasta", retmode="text"
) as handle, open(filepath, "w") as file:
    file.write(handle.read())

For batch operations, it would be better to acquire a NCBI API key and combine GNU parallel and Entrez Direct .
If you already have the accession numbers of ITS1 and ITS2 genes, put the accession numbers in a text file like "accession-number.txt". Otherwise search the ITS genes you want to download and then get the accession number:
esearch -db nuccore -query 'you query term[Gene Name]' |
  efetch -format acc >accession-number.txt

To download all the sequences in fasta format.:
parallel -a accession-number.txt -j8 'efetch -db nuccore -id {} -format fasta >{}.fa'

Here -j8 is used to control the number of jobs run in parallel, it should be less than the limit of 10.
To get related taxonomy (for eg, taxonomy rank, scientific name, taxonomy id, genus and species) of each sequence:
body='BEGIN { OFS = "\t" } { print "{}", $0}'
paralle -a accession-number.txt -j8 "elink -db nuccore -id {} -target taxonomy -name nuccore_taxonomy |
  esummary |
  xtract -pattern DocumentSummary -element Rank ScientificName TaxId Genus Species |
  awk '$body'"

Thank you! I thought, that they do not allow to make more than 10 requests per second for an API key regardless the fact it is parallel or serial. Which is real restrictions?

I think the real restriction is 10 requests per second. You can try -j with a number larger than 10 to see if there are any errors.

To get taxonomy info without awk:  paralle -a accession-number.txt -j8 "elink -db nuccore -id {} -target taxonomy -name nuccore_taxonomy |   esummary |   xtract -pattern DocumentSummary -pfx '{}\t' -element Rank ScientificName TaxId Genus Species"


text:
Multifasta protein alignment
I have around 4000 orthologous proteins for each of 5 species, I'd like to align them and to concatenate alignments. Which tools can I use? And how?
Marco

4000 isn't much. You'll need to start with an aligner. For aligners the common options are in no particular order:

Clustal Omega
MAFFT
FAMSA
Muscle v5

The current one that 'outperforms' the others strictly on a combination of two published benchmarks comprising 'accuracy' and 'speed' is MAFFT. FAMSA isn't relevant here unless you've a really small computer, like an old laptop, it excels at speed.
However, the situation could have changed recently with the  BioRXIV release of Muscle v5. Please note, I have not used it. Here's the Abstract:

On Balifam, ensembles generated by MUSCLE are shown to align an
average of 59% of columns correctly, 13% better than Clustal-omega
(52% correct) and 26% better than MAFFT (47% correct). The ensemble
bootstrap is applied to a previously published tree of RNA viruses ...

On this benchmark they are suggest Clustal Omega is better than MAFFT. Its not been the run of play that is certain MAFFT rose to dominance in recent history. Des Higgins Clustal W and X was historically the aligner of choice by a long shot, then Muscle, then MAFFT ... now its context specific.
Muscle v5 Its difficult to really get to grips with a BioRXIV ms as an algorithm unless you're a specialist in aligner construction. Please note, the 'RNA viruses' is complicated. As an algorithm this author has kudos and has extended that kudos beyond Muscle. The BioRXIV would have been subject to peer-review, but I am not aware Muscle v5 has been accepted for publication.
Historically, Muscle made early advances in accuracy which were duplicated by Clustal Omega at a much later date (2017?).
Essentially, here your choice is between MAFFT and Muscle v5 for this data set (context specific).
concatenate this is less clear but it appears you are performing 4000 alignments with 5 taxa in each? You'll need a data pipeline.
Trimming and concatenation are easily performed in BioPython using BioAlignIO but you would need coding because that can't be done manually. Computationally, BioAlignIO objects are annoyingly expensive - they are SeqIO objectives inside a MultipleAlignment object. However in your case it matters not a jot and would the recommended route all things considered (unless your coding expertise is Java).


text:
How to make a network overlap of gene lists across 3 conditions like this image?
From 3 gene lists (lung, nasal, gut), how to make this gene overlap network plot?
This image is figure 4a from: https://www.nature.com/articles/s41591-020-01227-z

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

I'm the author of gravis, an open source package for interactive graph visualization in Python.

It can visualize directed and undirected graphs (=having edges with or without arrows) as well as multigraphs (=having more than one edge between a pair of nodes).
It allows to control visual properties such as node sizes, colors, shapes or edge widths and color with data.
It recognizes graph objects from various network analysis packages such as NetworkX, igraph, graph-tool and others.
The output is based on web technologies (HTML/JS/CSS), so it can be viewed in a browser or Jupyter notebook. There is also support for exporting static images (JPG, PNG, SVG).
It comes with many examples on the documentation site.

Here's a quick attempt to reproduce a part of the figure you've provided:
import gravis as gv
import networkx as nx

s1 = 40
c1 = 'white'

s2 = 20
c2 = '#ffd703'

s3 = 20
c3 = 'orange'

s4 = 20
c4 = 'red'

g = nx.Graph()
g.add_node('Lung', color=c1, size=s1, border_size=1)
g.add_node('Nasal', color=c1, size=s1, border_size=1)
g.add_node('Gut', color=c1, size=s1, border_size=1)

g.add_node('MDFI', color=c2, size=s2)
g.add_node('LOXL1', color=c2, size=s2)
g.add_node('CNBP', color=c2, size=s2)

g.add_node('C3', color=c3, size=s3)
g.add_node('CA2', color=c3, size=s3)
g.add_node('SEPP1', color=c3, size=s3)

g.add_node('TFF3', color=c4, size=s4)
g.add_node('PTMA', color=c4, size=s4)

g.add_edge('Lung', 'MDFI')
g.add_edge('Nasal', 'LOXL1')
g.add_edge('Gut', 'CNBP')

g.add_edge('Lung', 'C3')
g.add_edge('Nasal', 'C3')

g.add_edge('Lung', 'CA2')
g.add_edge('Gut', 'CA2')

g.add_edge('Nasal', 'SEPP1')
g.add_edge('Gut', 'SEPP1')

g.add_edge('Lung', 'TFF3')
g.add_edge('Nasal', 'TFF3')
g.add_edge('Gut', 'TFF3')

g.add_edge('Lung', 'PTMA')
g.add_edge('Nasal', 'PTMA')
g.add_edge('Gut', 'PTMA')

gv.d3(g)

Resulting output in a web browser:

I've moved the nodes manually a bit so that the layout resembles the provided image. Once everything looks as intended, a static image can be exported with a single button click on the right-hand side.
There are many other ways to produce static graph visualizations, e.g. programmatically with other Python libraries mentioned before, or interactively in standalone tools such as Gephi, Cytoscape or Tulip.

You could use OmniGraffle to make this figure.
Objects in OmniGraffle have so-called magnets to which you can attach lines, which can join two objects.
Once you have objects joined, you could group them and position them as desired.
As you drag an object, any attached lines will also change their position as needed. You can drag groupings of objects (e.g., conditions or other categories of objects), which does the same thing.
You can export an illustration to PDF or SVG format from OmniGraffle, which will let you do any further work in Adobe Illustrator or other vector illustration tools, or share the illustration with colleagues over the web, for instance.

This question is a bit difficult to fathom, but essentially I think the answer you are look for is called network theory. Its quite simple, but alas not fashionable. You'll find an introduction here and try and avoid 'directed graphs' because that gets complicated.
The alternative way is to use Google fusion tables, which incorporates network analysis, which is really easy to use and investigate. Its very interactive. You'll find more information here,
https://support.google.com/fusiontables/answer/2566732?hl=en
Thats it and nothing more, undirected network theory is pretty simple. Its does not tell you very much (in my opinion), but it makes pretty diagrams.
Please note there is an analysis called 'Graph theory' and is not needed here. Whilst this has some similarity to 'network theory' it is distinct and the two shouldn't be confused.


text:
Issue with iCTNet2 installation
Has anyone encountered this same installation issue of iCTNet2?
I encountered this issue for installation of iCTNet2 via Cytoscape App Manager and similarly via installation from website download/install page.
I tried both installations via Cytoscape and via website takes a long time, when both is still installing even 4 hours later, I cancelled the installation. Does it really take that long? Or I should just wait for more than 4hrs to install this?
Any recommendations please? Thank you.

Answering this (2 years later) in case others stumble upon it. Apparently, the authors of this app are no longer supporting the server-side components. I was able to install the app, but not connect to their servers nor run the app. And the documentation link to their website is a dead link now as well.
So, we've removed it from the app store.


text:
How to identify most important Amino Acid residues in a given protein?
Continuing from the title, in my context, important means that changing the set of AAs (either or one more many) changes the domain and/or fold and/or function and/or family etc.
I  noticed that there are two types of methods Protbert (protein language models, sequence-based methods) and/or Alphafold (structure (tertiary) based methods); one needs both methods to "reallY" identify the essential set of AAs in a given protein?
Thanks

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.

Sorry. This question is ill posed as it completely depends on context —namely is simply a question of what are the main roles AA residues can take. Core residues when mutated may cause and increase in Gibbs free energy of folding (i.e. less energy is released when folded from a magically linear protein form), these are key to Structure as seen in PDB and predicted by Alphafold. Catalytic residues are key to activity. Channel gate residues are key to permeability. etc etc etc etc

Do you have many examples of your protein from different species? Can you look at alignments of the various homologs?


text:
Could someone please help me understand the outputs for HapFlk
Could someone please explain the output file for HapFlk?
My input file was a ped file which has three different populations. Two with same number of snps by imputation and third one is an outgroup population which has the same number of snps but some of them are marked as 0 since it's missing and unimputed.
I ran HapFlk with only the --file and --outgroup option so the output would be Flk not the LD based HapFlk and outgroup would be leftout after building the tree
When I checked the flk result after running,
the columns were rs, chr, pos, pzero, flk, pvalue which I could understand but some of the positions have

nan for pvalue
nan for flk and pvalue
1.0 for pzero and nan for flk and pvalue
0.5 for pzero and 0.0 for flk and 1.0 for pvalue

I'm not sure how to interpret these positions cause the positions aren't missing in input ped file.
Please help me

crossposted https://www.biostars.org/p/9526800/

It’s not against the rules, but it’s typically seen as polite to comment that you’ve posted elsewhere, to avoid people duplicating answering effort.

Regardless, you are likelier to get a good response to this at Bioinformatics SE :).

I'll keep that in mind thanks


text:
Elegant Way to Extract Gene Symbols from FASTA Headers
I've downloaded the RefSeq transcript sequences (FASTA format) for the GRC38 patch 14 build of the human genome and would now like to extract the gene symbols (only) in order from the headers. I'm confident I can do this using grep and regular expressions--the trick will be coming up with an appropriate pattern. (The gene symbol appears to be surrounded by the last pair of parentheses in a header, but I'm not sure this is absolutely true.) However, I'm wondering if there is a more elegant way to do this--e.g., via the NCBI EDirect utilities or SeqKit--that I'm overlooking. (I have a list of the accession.version numbers.) Any thoughts on this would be greatly appreciated.

It would be helpful to see a sample of the file or a direct link to the file in question. Is there a reason that you want to do this from the FASTA header (which is probably not so strictly curated) rather than from a more traditional annotation file format such as GFF3?

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.

We can't help you parse data if you don't show us what you need to parse. Please [edit] your question and show us a few lines (preferably just fasta header lines) of your file and the output you would want from that example.

We can't be sure without seeing your file, but if what you really want is the text between the last parentheses on lines that start with >, then you can do this with GNU grep (the default on Linux) with:
grep -oP '^>.*\(\K[^)]+' file.fa

The -o means "print only the matching portion of the line" and the -P activates Perl Compatible Regular Expressions (PCREs). The regex looks for a > at the beginning of the line (^>), then as many characters as possible until the last parenthesis (.*\(), and then forgets everything matched up to this point with the \K. This means that the -o will only see things after the \K, a PCRE feature. So finally, we look for the longest stretch of one or more non-) characters.
If you don't have GNU grep, you can do the same basic thing with perl or sed:
sed -n 's/^>.*(\([^)]*\).*/\1/p' file.fa
perl -lne '/^>.*\(([^)]+)/ && print $1' file.fa

Super helpful. Thank you!

You're very welcome, @MarkPauley. If one of the answers here solved your issue, please take a moment and accept it by clicking on the checkmark on the left. That is the best way to express your thanks on the Stack Exchange sites.

There are probably better ways of doing this with Python or another tool, but here's a quick method that works in bash (tested on Linux) assuming that the Gene Symbols are contained in parentheses, i.e. "(Gene Symbol)":
pat='\(([^)]+)\)'
while read line; do
    if [[ "$line" =~ $pat ]]; then
        echo ${BASH_REMATCH[1]};
    fi;
done < human.1.rna.fna

Thanks. Unfortunately, there can be more than one parenthesized expression in the header, so this pattern is too permissive. In every record I've looked at, the gene symbol is in the last parenthesized expression, but I don't know if that is guaranteed to be true.

Could you maybe use another file provided by the RefSeq ftp such as this one? It may have more genes than in your FASTA file, I'm assuming not fewer, but you could always do a quick check using the regex output as a comparison by sorting both files, and performing a diff command in Linux.

You really, really don't want to do this in the shell. It is going to be orders of magnitude slower than anything else and is very hard to get right. For example, your entire loop can be done faster and simpler with grep '\([^)]*\) human.1.rna.fna. The loop only makes it slower and more fragile.

Relevant: Why is using a shell loop to process text considered bad practice?


text:
How to identify a low proportion cell subpopulation in the single-cell RNA-seq data?
The cell subpopulation that I am interested in only accounts for around 1.2% of the total cells. I have previous FACS experiments that sort out the subpopulation from the samples (using markers CD166+ CD49fhi CD104- Lin-, see images). May I ask if there are any ways to identify this cell population from single-cell RNA-seq data of the tumor epithelial samples using the marker information?

Run the usual steps including clustering and visualization via something like UMAP and then color the UMAP by these four markers. That assumes that the surface marker separation you see in flow holds true on a transcriptional level which is not necessarily the case. For example, hematopoietic surface-defined populations such as the Granulocyte-Macrophage Progenitor (GMP) are well-separated in flow but on transcriptional single-cell cell level span many clusters.
Alternatively, check the tool UCell which can calculate a score per cell based on a combination of positive and negative markers. You can then color the UMAP by this to see whether a single cluster stands out or whether this celltype contributes to many clusters ir is not even present at all. Include some marker combinations that you know are oresent and some nonsense combination that must be absent to have positive and negative controls. That will give you a feel whether that score for your cells is reliable.
This all assumes that the cells are even in the dataset. Depending on cell size and frequency (underrepresentation) they might even have dropped out during library prep.


text:
How to investigate the correlation between transcript and protein abundance?
Our study shows that gene A knockout will down-regulate cell cycle-related genes in some cell lines. My aim is to understand whether cell cycle-related genes are down regulated on the transcriptional or translational level upon gene A knockout.
I plan to compare the transcript and protein abundance of cell cycle-related genes between two groups of cell lines (1st group: sensitive to gene A knockout; 2nd group: not sensitive to the knockout). May I ask how can I do this? Are there any tools, databases and papers related to my goal? Thank you very much!

This sounds more like a wet-lab question than a bioinformatics one. You are basically asking about how to run quantitative PCR or simlar techniques to measure transcript abundance and whatever would be used (sorry, no idea) to measure protein abundance, right? Once you have your data, we can help process them, but we can't help on the experimental side.

There are a large number of issues with this question. The first issue with the question is we don't know what organism is being studied, I might guess at yeast (traditional cycle model) - but that would be weird because there's loads of studies like that on yeast. It is very unlikely to be yeast because the RNA vs protein will be understood in detail. The second is what knock-down is it? ... and so on. Then there is the unfocused issue of the question about a comparison between transcriptomics vs. proteomics and how to analyse the relationship between them.

To simply answer the question, papers - yeast. Method of comparison supervised machine learning.


text:
Inconsistent replicate numbers in RNA-seq
I have 3 groups for the RNA-seq analysis (Control, treatment A and treatment B). There are 2 replicates for control and treatment A and 3 replicates for treatment B (lost 2 replicates due to a mistake in sample handling). May I ask how will the inconsistent replicate numbers impact the results (e.g. DEGs)?

It will not do any harm other than a reduction in statistical power compared to more samples per group. Common analysis tools such as DESeq2, edgeR and limma will work on unequal sample size just fine.


text:
Different results of spearman correlation between TPM and FPKM
TPM and FPKM of RNA-Seq data form GDC TCGA calculated based STAR were retrieved, respectively. The correlation between a specific gene, e.g. HIF1A, and other genes were calculated based on TPM and FPKM, respectively. And the significant genes were filtered based on P value < 0.05. Here is the intersection of these two significant gene set shown in venn diagram. IMO, the two gene sets should not be obiviously different. Ang advice for this discrepance would be appreciated!

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.

This shouldn't be surprising that you see different correlations between gene expression data when expressed in different units.
To see why, let's look at how these units are defined.
Let's denote the absolute number of fragments we sequence from an RNA-seq measurement as $X_{gi}$ for gene $g$ in sample $i$, and $l_g$ as the length of gene $g$ in kilobases.
Then our different expression values will be:
$$
f_{gi} = \frac{ X_{gi} \cdot 10^6 }{ l_g \cdot \sum_h X_{hi} } \text{ in FPKM}
$$
or
$$
t_{gi} = \frac{ X_{gi} \cdot 10^6 }{ l_g \cdot \sum_h \frac{ X_{hi} }{ l_h } } \text{ in TPM}
$$
You convert from FPKM to TPM units like so:
$$
t_{gi} = 10^6 \cdot \frac{ f_{gi} }{ \sum_h f_{hi} }
$$
As you can see, within any single sample, $i$, this last denominator is fixed.
But it is not fixed across samples (i.e. $\sum_h f_{hi} \ne \sum_h f_{hj}$ if $i \ne j$).
This means that the ordering for a single gene's expression across samples is not preserved when you switch between units (i.e. if $f_{gi} > f_{gj}$, it is not guaranteed that $t_{gi} > t_{gj}$).
Because Spearman correlation is based on the order of each gene across samples in your dataset, and the order of genes in FPKM will be somewhat different from the order of genes in TPM, you will find different correlations between HIF1A and other genes when using the different units.
Here's a simplified, concrete example.
Imagine that your transcriptome is made of up 3 genes (A, B, and C) with lengths 1 kb, 2 kb, and 3 kb, respectively.
And imagine you're doing RNA-seq in 4 samples.
If your absolute count matrix is this
$$
X = 10^6 \cdot \begin{bmatrix}
  1 & 2 & 9 & 7 \\
  2 & 2 & 6 & 5 \\
  3 & 2 & 2 & 5
\end{bmatrix}
$$
Then the expression matrix in FPKM is
$$
F = \begin{bmatrix}
  166\,666.7 & 333\,333.3 & 529\,411.76 & 411\,764 \\
  166\,666.7 & 166\,666.7 & 176\,470.59 & 147\,058.82 \\
  166\,666.7 & 111\,111.1 & 39\,215.69 & 98\,039.22
\end{bmatrix}
$$
and in TPM its
$$
T = \begin{bmatrix}
  333\,333.3 & 545\,454.5 & 710\,526.32 & 626\,865.7 \\
  333\,333.3 & 272\,727.3 & 236\,842.11 & 223\,880.6 \\
  333\,333.3 & 181\,818.2 & 52\,631.58 & 149\,253.7
\end{bmatrix}
$$
If you calculate the Spearman correlation for genes B and C against A, then you get these values.

Gene
Correlation with A (FPKM)
Correlation with A (TPM)

B
0.316
-0.8

C
-1
-1

You can see that gene C is anti-correlated with A, but B actually switches signs and goes from moderately correlated in FPKM to strongly anti-correlated in TPM.
So your observation that the same comparison on the different units gives pretty different results is expected.
If you made a scatter plot with Spearman's $\rho$ from TPM along the x-axis and Spearman's $\rho$ from FPKM along the y-axis for all the genes that had a correlation p-value < 0.05, you'll probably have a lot of genes that also flip signs like I showed with this toy example.
This paper by Wagner et al. is probably a good paper to review. It can give some context about RPKM and how to think about count-based sequencing data. Correlation-based analyses on data like this are typically used as a first-pass approach to get an idea about how things are behaving. But it's not very precise, as you've experienced. If you're doing exploratory work, I'd be inclined to pick TPM units, look at the results you have, then try to come up with a differential analysis question that can be handled in a more statistically robust way.

Got that! Thanks so much!


text:
Is the potential skewing of FDR and ranking of gene set analysis with gene sets with hierarchical structure a concern?
When performing GSEA and many other gene set analysis, when gene sets with hierarchy such as KEGG collection or Reactome collection, the FDR could be skewed sometimes, and the ranking of gene sets by any metrics could be confusing because multiple subsets belonging to the same set would have similar ranking, pushing other gene sets down the list. For example, in reactome axon guidance (R-HSA-422475.6) contains semaphorin interaction (R-HSA-373755.1), which contains Sema4d in semaphorin signaling (R-HSA-400685.2). The enrichment of Sema4d gene set usually (though not necessarily) entails the enrichment of the other gene sets above it. That might lead to an over-adjustment of p-value if I am not mistaken.
I am just wondering if there are ways to mitigate it, or it is not a major concern that I should have?

This is actually a known problem in the correction of nominal p-values in over-representation analyses that use the Gene Ontology categories (commonly but misleadingly called "GO enrichment" in many publications).
Most currently used methods for p-value adjustment (like BH) assume that p-values are independent, which is not true for a case like the one you describe - a set contained by another set is basically by definition not independent. Even if you didn't have a hierarchical structure, you would still have many overlaps across all your genesets, which would still mean that you cannot assume independence. So in both cases, but I think especially in your case, a BH/FWER correction would be statistically wrong (and possibly overly conservative).
There are a few ways to avoid this: one way would be to prune the directed acyclic graph of your genesets/categories, retaining only the categories that are enriched and have no enriched children.
In your case, for instance, if "axon guidance" and  "semaphorin interaction" are significant but "Sema4d in semaphorin signaling" is not, you would stop at "semaphorin signaling" removing the result for "axon guidance".
This is kind of a strong correction, which is guaranteed to lower the multiple testing correction burden (unless you have no significant genesets).
There are however more refined ways to address this:
One of them is to remove from set A (parent) any genes that belong to a significantly over-represented/enriched set A' (child). This method is termed the elim method and was described in this publication.
The same publication discusses another method called weight which assigns a weight to a significant category based on the significance of its children.
Another method is implemented in the gProfiler tool and introduced in this 2007 publication, which is the g:SCS method. Namely it calculates p-value thresholds based on empirical, randomized p-values for sets of different sizes. According to the authors it is more conservative than BH but less than Bonferroni.
The only issue with these approaches is that there are, to my knowledge, no ways to apply them to your data independently, i.e. they are coded within the respective tools and do not exist as stand-alone packages for multiple test correction. I hope I'm wrong and someone developed these versions, but for the time being the only thing you can do is to try to implement your own versions of these methods.
Another relevant, more recent publication is the SetRank one, which describes a method to do GSEA with a p-value calculation that accounts for overlaps, and has its own package on CRAN. If you prefer to apply a correction to your own results (like the ones you get from fgsea), however, I do not know whether this package has methods you can use or recycle.

Great answer! Is there a recommended standard package in R to trim DAGs without that one codes that up yourself?

@ATpoint IMHO it depends on the type of test. For Over-Representation Analysis (ORA), which is the one that normally uses GO genesets, packages such as topGO (Bioconductor) and gProfiler2 (CRAN) already implement trimming. For GSEA - what the user was asking about - SetRank may be one choice, but I have not used it personally. I am partial to fgsea for speed, so I did code my own pruning functions some time ago. I believe they are a little crude as they remove entire categories based on the DAG, but don't really address overlaps.

Thank you! Both answers provide some very good insight in gene set analysis. Trimming based on DAGs seems very straightforward and easy to do given database such as Reactome provides them in various format already. While searching, I also found a package called GeneSetCluster which takes GSA output and correct the results for overlapping genes between gene sets. I haven't tried it yet, but this could be an addition the DAG trimming and SetRank.

Statisticians might kill me for it but I've always filtered long and redundant geneset lists to only contain what remotely makes biological sense in my setup. If I am investigating leukocyte transcriptomes I do not care about terms related to filtration processes in the kidney or signal transduction in the hypothalamus. I think the redundancy of common pathway databases in terms of genes contributing to many terms justifiers to eliminate certain terms and choose those that make sense and are not overly large (for example not more than 500 genes). Top-of-hierarchy terms like 'metabolism' eith more than 1000 genes are anyway unspecific and hard to interpret. That then reduces the multiple testing burden quite a lot.


text:
postgreSQL error: "ERROR: invalid input syntax for type numeric: NULL"
I'm new to SQL. I have a table with a column of numeric  type that include NULL values. I built the table structure in postgresql, when I try to import my table I get the following error:
ERROR: invalid input syntax for type numeric: NULL
I tried both
UPDATE MY_TABLE  SET MY_COLUMN = CAST(NULL AS NUMERIC)
and setting manually the default value of the column as NULL, but none of them are working.
Is there another way to handle this problem?

Hi and welcome to the site! We only deal with questions about bionformatics here and this looks like a very standard SQL question with no connection to bioinformatics so I am afraid it is off topic. That said, I suspect you just want UPDATE MY_TABLE  SET MY_COLUMN = NULL.

I’m voting to close this question because it is about SQL and not related to bioinformatics.

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.


text:
VFDB search output
I want to check virulence factor of bacteria protein sequence by using Virulence Factor Database and it shows No hit in VFDB BLAST search output. Will someone please explain what is meant by it?

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.

Is your question about why blast doesn't get a hit, or about why your protein isn't in VFDB? BLAST and VFDB both have documentation about how they work, can you try to investigate there?

if VFDB result show No Hit for a protein, so is it means it is non-virulent?


text:
Converting .bw files to .fq (fastq)
I am able to convert .bw files to .fq(fastq) manually in Shell, but I would like to automate the process coz I have hundreds of .bw file that I need to convert. Till now, I could think of following script, but unable to solve the problem. The .bw files are in H3K27ac folder.
import subprocess

path_bw = "~/Documents/H3K27ac"
import os

files = os.path.expanduser(path_bw)
for file in files:
    if file.endswith('.gz'):
        files.remove(file)
files = [x[:-3] if x.endswith('.bw') else x for x in files]

# for .bw files

for file in files:
    bigwigtowig = "bigWigToWig " + path_bw + file + ".bw " + path_bw + file +".wig"
    subprocess.call(bigwigtowig, shell=True)

    wigtobed = "wig2bed < " + path_bw + file + ".wig > " + path_bw + file + ".bed"
    subprocess.call(wigtobed, shell=True)

    getfasta = "bedtools getfasta -fi " + path_bw + "hg19.fa -bed " + path_bw + file + ".bed -fo " + path_bw + file + ".fa"
    subprocess.call(getfasta, shell=True)

    perl = "perl " + path_bw + "fasta_to_fastq.pl" + path_bw + file + ".fa > " + path_bw + file + ".fq"
    #os.system(perl)
    subprocess.run(perl, shell=True)

Manually, I used the following to convert the .bw files:

# bigWigToWig file1.bw fileout.wig
# wig2bed <fileout.wig> fileout.bed
# bedtools getfasta -fi hg19.fa -bed fileout.bed -fo fileout.fa
# perl fasta_to_fastq.pl fileout.fa > file1.fq

All the above packages were preinstalled via shell and works perfectly. hg19 is the genom that i am using. fasta_to_fastq.pl is openly available in github and helps to convert .fa to .fq.
MAin thing i require is to convert all the files from .bw to .fq and this is the only way in my research that i could find. I need all the names to be the same as output.  For example file1.bw shoud be file1.fq , file2.bw --> file2.fq and so on.
Do you have any ideas if this is possible or are there easier ways to convert it? It would be awesome to know your view point.

Can you please explain why you want to do this? It's a very unusual situation, and it's likely that there are better ways to approach what you want to do. It is often better to ask questions about your actual goal, rather than intermediate steps (see here for more information).

What you are doing makes no sense at all. A bigwig file contains intervals and then indicates the coverage of that interval. That coverage again comes from the alignments covering that region. There is no direct connection between the interval and the fastq file other than that the fastq file is the basis for the coverage calculation. Converting the interval to a fasta and then forcing it into a fastq file is technically possible, but has absolutely no practical meaning. You cannot recover sequencing reads from a bigwig, that information is lost once you convert your alignments (usually a BAM file) into a bigwig. You should open a new question describing the problem that this is actually based on.

I'll add slightly to this excellent answer and mention that if this is public data (highly likely), there is almost certainly a link somewhere to download the underlying fastq files.

Excellent answer that misuses the word "convert" like most people on these forums do. One does not convert BAM to (big)wig, one extracts/computes wig information from a BAM file. "Convert" kind of implies information content preservation (convert sam to bam, etc.), which as you say in your answer is not true when getting WIG from BAM.

The code is verbose and I ain't at all sure some of it will work, at least in two places.
The first point is os.path.expanduser ... admittedly I prefer pathlib, however that is definitely not how I understand that command. My understanding of expanduser is if ~ is being used. You definitely want to look at globb-ing for the opening part of the code.
from pathlib import Path
import os

mypath = Path(os.path.expanduser('~/Documents/H3K27ac'))
files = mypath.glob('*.bw') # or use ('*' + '.bw')

I'm pretty certain that is what you want and its a lot easier than lines 1 to 10 in your code. The files list will then go straight into the loop.
Also, please check
files = [x[:-3] if x.endswith('.bw') else x for x in files]

Comprehensions are never logical in Python but convention is always to place the if statement after the list of interest.
Thus the following is much closer:
files = [x for x in files if not x.endswith('.bw') else x[:-3]]

Note not tested, so this may throw a bug, but its better. However, using glob bing you don't need this statement.
The other issue is the sub.process call ... Normally they are written:
os.chdir(path_bw)
file1 = 'file1.bw'
file2 = 'fileout.wig'
subprocess.call(["bigWigToWig",  file1, file2], shell=True)

The way you've written is:

Really hard to read
I don't think it would work because of the absence of the [], certainly if you are using popen it will not work

Also check the following variant:
subprocess.call(["bigWigToWig",  file1, file2], text=True)

I know its seems a bit weird, but it very likely work (don't ask why)
I agree subprocess.call is a good idea in a Python pipeline. Issuing a hard path within subprocess call I'm not a fan hence using os to get Python to move into the cwd is a better approach.
A further point is the way you've used 'file' throughout your pipeline line, hmmmm... okay, but its difficult to read. Its not generally how a pipeline is written, but just saying.
nextflow is an increasingly popular pipeline language and you may consider learning that rather than going down a Python route (Python is very good for data science), of course this particularly question could be very easily solved in bash (but for some reason it ain't trendy).


text:
What tool I can use to map short-reads sequences to reference genome and get specific mapped size
I have about 300 90-bp sequences which I would like to map to a reference genome to make it longer to a 300bp sequence, wherein the 90bp is at the middle.
Anyone knows what bioinformatics tool I can use for this? Thank you!
Actual scenario:
The SNP sequences are from a SNP panel data (90bp sequence, SNP is at the middle). We would like to extend the sequence to 300bp, retaining the SNP at the middle so we can make a sequence capture panel of the SNP.

You mean the target capture panel of SNP markers?

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.

@envs_h_gang_5 I'm actually not very familiar, but something like 'genotyping by sequencing' kind of panel.

If you have a SNP panel then you have the coordinates on the genome. You can use tools such as bedtools slop to expand the coordinates to the desired length and then extract those regions from the target genome using bedtools getfasta.

Thank you! That's very helpful. The dilemma, however, is we'd like to extract the 300bp sequences from a different genome, and not from the reference genome of the SNP panel, therefore the coordinates on the SNP panel will not be exactly on our genome.

Then you need liftover https://genome.ucsc.edu/cgi-bin/hgLiftOver
I don't know if the genomes you're working on have mapping files though.


text:
Find enriched processes in modules obtained by WGCNA
I am doing a gene enrichment of different clusters obtained by WGCNA, but there are some of them in which I don't get enriched processes, is there any way to find processes to these modules, for example to associate them to orthologues of other organisms. If there are different ways, could you send me links on how to do it?

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.

Hi María José. What kind of gene enrichment are you using? Which databases are you using to infer functionality? What organism are you analyzing? If it is human I doubt there are other organisms more annotated than that.

Well, not all modules will have significant enrichment in gene sets such as GO terms or KEGG pathways. In my experience, the smaller modules are often not significantly enriched. They could be the result of some technical variation, they could contain genes localized to a particular genomic region (do you also calculate enrichment in genomic regions?) or they could be a particular gene family. WGCNA modules often are enriched in known pathways and processes, but that does not mean that all will be.

thank you so much for your clarifications @Peter Langfelder , I liked the idea of calculating the enrichment by gene regions, which packages would you recommend for this.


text:
Indexing the reference genome my process is killed
I am trying to use pbmm2 to create a reference genome index.
I'm on a cluster server, I don't know if that's the reason why it has been killed.
pbmm2 index hg38.fa hg38.mmi

You're probably running out of memory. Please Google before asking, this is a Linux question and not specific to bioinformatics.

@RamRS I have tried with ulimit -m or -v ~64GB; I don't know how many GBs would be needed for a human genome, previously it was indicated unlimited.

How many gigs depends on the program. How much physical memory do you have access to though? Given that pbmm2 is essentially minimap2 and minimap2 was developed by Heng Li, I'd say it should be pretty efficient, but try using minimap2's indexing operation instead of pbmm2's.

Are you running this job on the head node? It has limited resources, and is only supposed to be used to communicate jobs to the compute nodes. Seek out your server administrator and see if they have a wiki or other instructions on how to run jobs on their cluster.

@Pallie assuming that the cluster is running SLURM the queues/partition looking at /etc/slurm/slurm.conf

That Killed message is what you see when a process is killed using a KILL signal:

This is the signal that the out of memory (OOM) killer, the kernel process whose job it is to stop processes if they are using too much memory. Instead of allowing the system to hang, the kernel starts killing things that are hogging the memory.
I have never used pbmm2 myself, but I checked the documentation and it doesn't appear to have any options for memory management. This absence suggests that it doesn't really expect memory to be an issue, in which case it is surprising that you got that error.
Could it be that you are launching this on the host node of your cluster and the cluster is designed to work with a queuing system that sends jobs to the full cluster instead? Or perhaps the system was under very heavy use and other users' jobs had taken all the RAM leaving none for you? Essentially, it looks like you will need to run this on a machine with more RAM.
Alternatively, on the off chance that the tool is just badly designed and doesn't manage memory well, since all you want to do is make an index, you could try using minimap2 in case that works better. You can also not create an index and let pbmm2 create it on the fly, assuming it works like minimap2. However, if your system doesn't have enough memory to make the index, I really doubt it will have enough to create the bam either, bringing us back to: use a machine with more RAM.

I got the same problem on an 8GB machine, where there was only about 4GB memory left available. I was able to build an index on a bigger machine with more memory.


text:
macs3 to detect high coverage regions in long-reads
I am using macs3 to detect regions with high number of reads, similar to macs3's goal - detect chip-seq peaks - as it has been used for long-reads data, I would like to ask how to determine if it is doing a proper call or detecting spuriously peaks. The peaks are analyzed through IGV visualization.

Can you please provide more context around your study, especially examples of outputs / peaks? It's hard to answer this without knowing what data you have.


text:
Verify cell identity
We want to analyze from an article the sequencing data (3 rd generation) of the U2OS line, and we would like to verify if it is that cell line. I was wondering how would be the general procedure, or if there is any tutorial for it.
I was thinking about performing a STR profiling with ExpansionHunter or Tandem Repeats Finder, for example of CSF1PO or D13S317. I was also wondering if I would need to perform a Variant call with GATK before.


text:
Generating contigs from illumina Paired End Reads (Metagenome)
I'm a wet-lab biologist desperately trying to get to grips with some Bioinformatics.
I have some shotgun Illumina paired-end data from a gut microbiome project I'm working on. I'm considering how to go about generating contigs from reads. I've been recommended to use SPAdes which I'm reading up on at the moment.
Originally I was thinking of using FLASh however I was told it might not be suitable for metagenomic data.
My question simply is, what makes certain programs suitable for merging paired-end reads from a genome but not a metagenome? I find myself a little spoilt for choice (eg. FLASh/SPAdes/BBmap... )  with bioinformatics tools I'm always concerned I'm picking the wrong one.
Should I ideally be running my data through a couple of tools and comparing them?

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.

Is the source material metagenomics or culture isolates?

What is the purpose of your project? You have mentioned, "merging paired-end reads", but it's unclear why this would be helpful. If this is for quantifying species or determining proportions, then you'd be better off with a metagenomic mapping tool (e.g. kraken2).

NB: "merging PE reads" (such as e.g. FLASH does) is not the same as "sequence assembly from PE reads" (what MetaSPAdes does). The latter generates contigs, the former generates (some) merged shotgun library inserts. I think you want sequence assembly, i.e. a tool like MetaSPAdes.
Don't be tricked by the loose language of "merging" that is sometimes applied to the sequence assembly problem. It is much more complicated than that.
If you want to learn more about sequence assembly, I suggest you read the sequence assembly wikipedia page.
I would probably recommend MegaHit over metaSPAdes, because it requires much less memory / runtime, and lower frequency of errors.
But generally speaking, both programs are building different flavors of de Bruijn Graphs from metagenomic reads to generate contigs/scaffolds. This is the accepted method by this point for Illumina metagenomic data, because it is the method that is capable of handling the amount of small pieces of data in a tractable fashion. Alternatives such as overlap-layout-consensus graphs or string graphs just don't work as well for the metagenomic problem.
You are probably looking for methods that talk about complicated, computationally complex assembly algorithms and graphical methods for summarizing lots of short reads. You are not looking for simple read processing/alignment tools like BBMap/BBTools or FLASH. To reiterate, simple read merging tools will not give you assembled contigs.

Thank you, this had helped a lot, I'm definitely struggling with some of the subtleties of data processing. My end goal is to do Taxonomic and Functional analysis in MEGAN, but I've been told that it's only designed to use R1 data. I thought I'd attempt to make the most of the sequencing data that I do have to get into a format suitable for MEGAN. Would using  MegaHit/metaSPAdes mean I wouldn't have to put my data through DIAMOND ?  I'm still very new to this and just making sure I have all the steps correct

I am not sure what you mean by only R1 data. I don't think that assembly will make DIAMOND work any better, though it may be helpful for some MEGAN workflows. For examples of both assembled/raw shotgun data for MEGAN+DIAMOND, I'd suggest this ref: https://currentprotocols.onlinelibrary.wiley.com/doi/10.1002/cpz1.59

Thank you, I'm still very much getting to grips with this.  I originally did 150 paired-end sequencing on my samples. I think DIAMOND/MEGAN only accepts a single input file (What I meant by only using Read1 data, I wasn't being clear there.) I thought trying to merge overlapping reads would give me longer DNA sequences for Taxonomic binning.

@dunc4n I'll first note as a tip for the future that having this intention stated in your question would make it a lot easier to understand what you need. I suppose that it's formally possible that read merging would improve things slightly, in practice I am not sure that it makes any difference. I would follow the example of the reference that I linked in my earlier comment, or really any other paper that uses DIAMOND and MEGAN, to understand what to do to preprocess data for those tools. Is there a reason why this won't work for you?

I'll keep that in mind in the future thanks. I don't see any reason why it would work. What had me confused was when I mentioned my plans I was advised that FLASh wouldn't be suitable to generate contigs from my metagenome samples. I had planned to generate contigs from Read1 + Read2  and potentially cat on the output of the unmerged read by FLASh or align them separately and join them post assembly.   I wonder if when I mentioned it originally They thought I was intending to rely on it for contig generation and assembly.   Thank you again for your input it's been really helpful !

@dunc4n  yeah, so I think that the source of the confusion is that merged reads aren't "contigs". Merged reads may represent "contiguous sequences" but the word "contig" has a rather specific meaning in genomics. I'm still not totally clear on the workflow here, but it seems like the issue is resolved at least.


text:
Downloading genomic protein files from accessions in Python
I am trying to download the _protein.faa.gz files for genomes given their accession numbers through Python. Ideally, I would like to do this without third party libraries. Essentially what I have is a list of just the GCA or GCF accessions. The issue with the ftp site is that it includes the project name along with the accession, which I do not have ahead of time. It would be perfect if I could run downloads from within Python that matched a pattern something like this:
https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/145/295/GCF_000145295.1_*/*_protein.faa.gz
Is it at all possible to do this? Or is it possible to run something similar via efetch in Biopython?
Thanks in advance!

Concept code only
from Bio import Entrez

Entrez.email = "m@M__"
handle = Entrez.efetch(db="protein", id="GCF_000145295", "GCF_000145294", "GCF_000145293", retmode="text",rettype="gb") # check the quotes
records = Entrez.parse(handle)
handle.close()
for record in records:
    print (record)

You can write a loop, flatten it and then insert it to the id= tag.
mylist = ['GCF_000145295', 'GCF_000145294', 'GCF_000145293']
flattened = ','.join(mylist)

The above answer is right - here are some more links that describe what you want to do in more detail:

https://eeob-biodata.github.io/BCB546X-python/06-biopython/
https://medium.com/@insidedna/biopython-to-retrieving-fasta-sequences-from-ncbi-databases-20e35772a25

You're going to want to use the BioPython package as it has functions that do exactly what you're looking for.


text:
Different line length in fasta file
I am currently using VEP for variant annotation.
I am facing an error as below:
[E::fai_build_core] Different line length in sequence 'Pn9'

I understand there is an issue with the difference in line length of Pn9 in the  fasta file. However, the sequences of fasta file all have different line length. I don't get why there is an error specifically on Pn9.
I have tried using both snpEff and Annovar but doesn't work. Any thoughts on this would be really appreciated.
I am attaching the length of sequence here :

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.

This could be a bug in the fasta file reader component of fai_build_core. The fasta file format is very relaxed, and allows some odd constructions, for example:
>cb3c6a5f-17dd-4cd6-a9ca-02124d2a32de
AAGAAATTGACCTCCGTAACCCCAAAGGAGACAAAAAGAATGCTGGCTTTACTCCC
CAGGAGCGCCAGGTTTTGGGGAACTGCCTACAAGCTGTACCATTGGCTGACAGCTT
TTAGGGCATCTACCCCCAACACTGCTTACACTTACACTTTCTGGACTTACATG
ATGAATGCCCGCTCTAAGAATGTTGGTTGGCGCCTTGATTACTTTTTGCTTTCC
CACTCTCTTTTTACCTGCATTGTGACAGCAAGATCCGGTCAAGGCTCTTG
GCAGTACACTGTCCCATCACCGCTACCTGTGATACCCCTCCCTGCAGTAGCTT
CCTGCCTGGGAGATGGCTCTCTCTGCAGAAGTCTGGTGTTTTAGCC
TTCAGGTGTTTGGTTTGTATGTGCTCCCTCATTTAAACATTAAACCAAACTT
CTGGTTTCCTTTAGACAATCCAAGAGAAATAAAAGGCCCTAC
TTTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA

A few programs will complain about fasta files that look like this, because it usually means that data has been corrupted. These programs only accept fasta files that have sequences with identical line lengths for all but the last line:
>cb3c6a5f-17dd-4cd6-a9ca-02124d2a32de
AAGAAATTGACCTCCGTAACCCCAAAGGAGACAAAAAGAATGCTGGCTTTACTCCC
CAGGAGCGCCAGGTTTTGGGGAACTGCCTACAAGCTGTACCATTGGCTGACAGCTT
TTAGGGCATCTACCCCCAACACTGCTTACACTTACACTTTCTGGACTTACATGATG
AATGCCCGCTCTAAGAATGTTGGTTGGCGCCTTGATTACTTTTTGCTTTCCCACTC
TCTTTTTACCTGCATTGTGACAGCAAGATCCGGTCAAGGCTCTTGGCAGTACACTG
TCCCATCACCGCTACCTGTGATACCCCTCCCTGCAGTAGCTTCCTGCCTGGGAGAT
GGCTCTCTCTGCAGAAGTCTGGTGTTTTAGCCTTCAGGTGTTTGGTTTGTATGTGC
TCCCTCATTTAAACATTAAACCAAACTTCTGGTTTCCTTTAGACAATCCAAGAGAA
ATAAAAGGCCCTACTTTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA

One way to work around this input error problem is to get rid of all line breaks in the sequence of a fasta file, but then you run the risk of programs not properly being able to read all the characters in a line. A better approach is to feed the fasta file through a less restrictive data cleaning program and convert the lines back to all have a consistent length.
There's a good chance this is not the problem you are experiencing, but without more information about your input data, all we have to go on are the error messages and the sequence lengths [not line lengths] you have provided.

Thank you so much for your reply. I was able to rectify the error. I was doing some unwanted steps while zipping the fasta file.


text:
Get number of reads with a single, (almost) exact match to the full length of a reference sequence
I can't find an answer to this in previous questions, so hoping someone can help me now.
We have nanopore sequenced a PCR product, and I've filtered our reads to +/- 10 bp of the expected product length (~4,400bp), leaving ~55,000 reads.
We now want to work out how many of the reads correspond to each of 6 different potential arrangements. The section we're interested in is in the middle of the 4,400 bp PCR product, and each different potential arrangement is ~1,500 bp. We have reference sequences for each of the 1,500 bp arrangements. Basically, we want to know what portion of our starting cells had each of the different arrangements of genes in this operon (e.g arrangement 1=30%, arrangement 2=15%, etc.).
This seems like it should be a really easy thing to do using read mapping to count alignments for each reference separately, but I just can't get my head around it. No matter what options I've tried with minimap2 and samtools, every read maps every time on every different arrangement. I think the problem is that the genes in our reference sequences are all present in every read, just in different arrangements. So we need to only count alignments which align to the whole 1,500 bp reference sequence, not partial hits or hits with huge gaps in them.
I have tried the minimap2 options -r, -G and --secondary=no to restrict the output to single alignments with no large gaps, and combined that with samtools view -m 1000 to output only alignments longer than 1,000 bp (this is far stricter than it needs to be to only output alignments to any of the different potential arrangements). But still, every time, every read is mapping to every reference, usually multiple times.
What am I doing wrong?! Am I totally over-complicating this, or doing something obviously stupid? Almost feel like using blastn might be better at this point...
TL;DR: how do I use minimap2 and/or samtools (or any other tool better suited to this) to count how many of our ~4,400 bp reads contain only full-length, almost exact matches to our 1,500 bp reference sequence, excluding any partial hits?
My latest attempt (broken into individual steps, though I'd usually pipe the samtools steps together, also not bothered with samtools sort at this point):
minimap2 -ax map-ont -r10 -N5 --secondary=no --no-long-join reference.fasta reads.fastq > reads.minimap.sam

samtools view -m 1000 -b -o reads.minimap.bam reads.minimap.sam

samtools -c reads.minimap.bam

Because minimap2 is a minimiser-based mapper, it's less useful for determining mapping accuracy to single-base precision.
Where precision matters, you can use LAST, bearing in mind that it calculates local alignments, and can report multiple hits per query/reference pair. The usual approach I use for LAST mapping is as follows:

lastdb -uNEAR -R01 reference.fa reference.fa
last-train -Q 1 reference.fa <(cat reads.fastq | head -n 40000) > train.mat
lastal -p train.mat reference.fa reads.fastq > reads_vs_ref.maf

That final MAF file can be inspected for mapping. The default output produced from LAST is not very human-readable, so I've written a perl script that converts MAF results into a CSV format with reference-relative statistics that I have found useful:
maf2csv.pl reads_vs_ref.maf > reads_vs_ref.csv
This can be inspected in a spreadsheet viewer (e.g. LibreOffice Calc, or Excel) to look at the mapped locations, mapped fraction, identities, etc..
Alternatively, the MAF results can be converted to SAM/BAM format using maf-convert (bearing in mind the previously-mentioned caveats about multiple hits):
maf-convert sam reads_vs_ref.maf | samtools view -b -T reference.fa | samtools sort > reads_vs_ref.bam


text:
AttributeError [in fastStructure's chooseK.py]: module 'vars' has no attribute 'insum'
I've been working with the fastStructure program and am on the step of analyzing model complexity using the provided chooseK.py script. I have been running into the following error:
Traceback (most recent call last):
  File "chooseK.py", line 109, in <module>
    bestKs = parse_varQs(files)
  File "chooseK.py", line 54, in parse_varQs
    Q = Q/utils.insum(Q,[1])
AttributeError: module 'vars' has no attribute 'insum'

I have run into 'AttributeError' messages in the past and have attempted the solution of renaming files in my working directory that may be shadowing the real module (similar to: https://stackoverflow.com/questions/36530726/using-numpy-module-object-has-no-attribute-array), but this hasn't fixed anything. Reinstalling Numpy/Scipy/fastStructure has also been a dead end.
Any suggestions are much appreciated! Happy to provide more info as needed.
-A

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.

did you build the library extensions? The module vars.utils are imported as utils, so the attribute insum is undetected. https://stackoverflow.com/questions/46312470/difference-between-methods-and-attributes-in-python

Could you try loading fastStructure via conda?
conda create --name faststruc python=3.9 # there can be issues with 3.10
conda activate faststruc
conda install -c bioconda faststructure

It will then ask you whether you want to install numpy and scipy automatically. Furthermore all dependencies are installed and connectivity will have already been tested. You can easily delete everything afterwards by deleting the environment faststruc so keeps everything clean.

The code in question is here:
https://programtalk.com/python-examples/vars.utils.insum/
Its a numpy calculation, what they are doing is navigating Python by leveraging the power of numpy, which is very fast.
Code ...
bestKs = []
for file in files:
    handle = open(file,'r')
    Q = np.array([map(float,line.strip().split()) for line in handle])
    Q = Q/utils.insum(Q,[1])
    handle.close()
    N = Q.shape[0]
    C = np.cuemsum(np.sort(Q.sum(0))[::-1])
    bestKs.append(np.sum(C<N-1)+1)

return bestKs

The code looks good, they know their stuff.


text:
Establishing Quantitative Links Between Regulatory Activity Marker H3K27ac ChIP-seq and Differential Gene Expression, Without Hi-C or PLAC-seq Data
I am working on correlating H3K27ac ChIP-seq data, which indicates gene expression-related activity in regulatory regions, with RNA-seq gene expression data to establish quantitative connections between gene expression and regulatory activity marker H3K27ac. Specifically, I am comparing two distinct cell types, resulting in differential gene expression (DG = log2FC(G2/G1)) and differential H3K27ac (DH = log2FC(H2/H1)) based on ChIP-tag count differences. My approach involves the following steps:

Identify Active Differential Peaks: Select H3K27ac peaks with |DH| above a certain threshold.
Associate Peaks with Genes: For each active differential peak, identify the closest 2-3 genes based on genomic distance.
Intersect with Differentially Expressed Genes: Determine the overlap between these associated genes and those with |DG| above the threshold.

This method assumes that differentially active H3K27ac peaks primarily regulate nearby genes. However, I am not sure if this a good assumption, as active peaks might regulate genes in much farther proximity than nearest two or three genes. In the absence of any additional data such as Hi-C or PLAC-seq, which can explicitly link promoters and distal enhancers, are there alternative methods to narrow down a list of genes that maybe associated with specific active enhancer peak region? I would appreciate any critiques to this approach for connecting H3K27ac activity and gene expression.


text:
How to extract DNA sequence from browser track input files (BigWig, bed etc.) files
A BigWig file for something like Chip-seq, or ATAC data has the following format for each line i:
chrom_i - start_i - end_i - value_i

I need to convert for each line the (start_i -> end_i) reference to the actual DNA sequence 'ATATCCTAG...' so that at the end I have a list of DNA sequences like 'ATATCCTAG...' and the Chip-seq activity value associated with that sequence segment. What are the ways to do this? For example, what tool(s) can I use to generate the (sequence, value) pairs from browser tracks for the Chip-seq or ATAC etc. And in fact, going from BigWig browser track files for Chip-seq the most direct way to produce this kind of (seq,value) data? Thank you.

[Here is the link.][1]. Instructions are very clear at the website
[1]: https://bedtools.readthedocs.io/en/latest/content/tools/getfasta.html

Ah, perfect @Supertech! Thank  you, both answers helpful... the link included here is the exact answer I was looking for: bedtools getfasta converts index range start_index-end_index into it's equivalent sequence of nucleotides indicated in the input *.fa file. I just need to figure out how to create the "fa" file for an entire genome like mm10, where it's usually in chromosomal chunks (i.e. Chr2.fa).

I believe each chromosome is independently numbered. If you want to generate a single genome file with all the chromosomes, download individual chromosomes and simply concatenate them. For example using Linux command "cat". >cat chr1.fa chr2.fa ...

It looks like they may have single file. I have not downloaded to check it though. It has been long that I did not work with mouse genome. http://hgdownload.cse.ucsc.edu/goldenpath/mm10/bigZips/

Yes, I am trying to think why/how that would work: I guess the *.fa file has the chrX indicator at each line which would be matched along with the index (I assume a complete index needs to include chromosome number as well)! Thank you again.

Try samtools faidx. Check bedtools manual. https://bedtools.readthedocs.io/en/latest/content/general-usage.html

Bedtools can take a list of coordinates and a fasta and return a multifasta of the specified regions.


text:
Tool for cross referencing a newly found MOTIF to database of known MOTIFS
I need a tool or a function I can use in my code (R, or Python) that I can cross reference a MOTIF against known MOTIFS, a function that will take as input a MOTIF (a probability weight matrix, PWM which is a 4 by m matrix for an m-length motif) and spit out the closest -n- known motif(s). Is there such a function/tool or something close to it? Thanks much.

A popular tool for doing this is tomtom, part of the meme package.
You can try it via interactive web form here.
Once installed locally it can be run via command line
Additionally, it is available in the Bioconductor package memes as explained in Motif Comparison using TomTom

Thanks @malcook: TomTom seems like the tool.  The list of required libraries/tools to install first before tomtom can be installed is quite lengthy and even after that the installation is riddled with errors in the "make" step for me. Any hints on installation on Linux/Ubuntu?

Hi @ZebraFish - I have no hints - you are right: the meme suite is a bit of bear to install.
The R/BioConductor package 'universalmotif' also addresses this purpose and is  easier to install.  I have not (yet?) used it for such, but have been recently impressed by its documentation specifically that on "Motif comparisons and P-values".

Thanks @malcook: After much haggling and hacking around I discovered the simple way to install TomTom is actually through "conda install -c bioconda meme". This puts a slew of meme command line executables into a conda env. directory including TomTom. This simple method of installation was unnecessarily hard to find  if you are new to this field. Just looking at the meme suite webpage It's a curiosity why such simple tools (from a computational person's perspective motif search against a database of motifs is a very simple task) can be so obfuscated.

Answer from @m, converted from comment:
Hi @ZebraFish Motifs package Biopython
https://biopython.org/docs/1.76/api/Bio.motifs.html


text:
Can FASTA data be submitted directly into eCTD submission to the FDA?
We are getting sequence data from microbiota samples as part of a FDA submission identifying changes in Micobiota from skin samples over time and compared with different treatments.  the data will be first presented in FASTA formate then analyzed for OTUs and differences.  The submission has to follow eCTD format.

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.


text:
Evolutionary tree from SNP data
I would like to make an evolutionary tree using the SNP data from DArTseq. Are there any software or packages available for this purpose? I came across the software SNphylo, but its home page is not updated. Your suggestions would be appreciated.

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.

Welcome @Lama You would need to explain DArTseq and the data set size

At present I would simply use iqtree
iqtree -s file -m GTR+I+G -T 6 -B 1000

This will parallelise across 6 CPUs and results in 1000 bootstraps, whilst GTR is the mutation model, which models invariant sites under a 4 category gamma distribution. This is a robust model.
You can view the tree in FigTree.
The DNA file is in fasta format
> seq1
AGACGT
>seq2
AGAAAT

If there's an 'unknown' format bug let me know because I only work in  data pipeline so I forget if I've shifted the format.
The result will be in the treefile and is a professional quality tree, complete with bootstraps.
The I parameter here could omitted for this dataset because I assume all none of the sites are invariant.

@M_ Thank you for you answer. DArTseq is one of the cheap and easy but efficient genotyping-by-sequencing platforms which allow genome-wide marker discovery through restriction enzyme-mediated genome complexity reduction and sequencing of the restriction fragments (Edet et.al. 2018). I am attempting to construct a phylogenetic tree using the results of DArTseq rather than the FASTA sequences.

@Lama its a simple format shift. Maybe conversion of DArTseq to fasta is a good question to raise? At the end of the day the theory is universal if the output is nucleotide data. We'd need to see the new type first, but it's likely very easy. You can build trees from anything and everything - its the underlying theory that counts.


text:
How can I delete these lines in my fastq file?
Hi guys =) I'm sorry if this is a repeat but I haven't been able to find an answer using the search field or google.
I am trying to edit a fastq file that has a corrupt line. The read in question is VH00362:5:AAAV5KTHV:1:1210:64850:6869 where the read is just a massive repeat – not the expected 110bp. I am trying to remove two lines from this file and asking around  I've been told to use sed as in
gzip -c file.fastq.gz | sed "45500973,45500974d" > new.fastq.gz
but when I run wc -l  to check the number of lines the results are incoherent, as if it has deleted much more than the two lines I asked. Another issue that I'm seeing with this file is that every time I run gzip I get a invalid compressed data--format violated error, which tells me the file is corrupted but I'm doing my best to salvage whatever I can from the file as it is an irreplaceable sample.
Is there another way to simply delete these two lines? Any tool or command you could recommend? Any help is deeply appreciated.

Does gzip -k file.fastq.gz report an error? If so, I would first try to gunzip the file, or zcat file.fastq.gz > file.fastq and then operate on the decompressed one. I am not sure, but I suspect that a corrupted archive might not cat as expected and I would feel safer getting as much of the file decompressed as I can first. I am not even sure you will get anything, please try and report back.

This is not the thing I do these days but ...
gunzip -c file.fastq.gz | sed "45500973,45500974d" > new.fastq

You were using gzip instead of gunzip, and that the source of the error. Finally the out file isn't compressed, so you probably want:
gunzip -c file.fastq.gz | sed "45500973,45500974d" | gzip > new.fastq.gz

That should work (I tested it and yes it works).

A fastq file has four line-separated fields per sequence. If you suspect that a sequence is corrupted, then you would want to delete all four lines related to that sequence. If you just delete the first two lines (sequence identifier and raw sequence letters), the file will most probably throw an error down the road. In your case, the sequence would look something like this -
@VH00362:5:AAAV5KTHV:1:1210:64850:6869
GATTTGGGGTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65

So, the code to delete these four lines from a fastq.gz file would be -
gunzip -c file.fastq.gz | sed "45500973,45500976d" | gzip > new.fastq.gz

This sed command specifies a range of lines and says that delete lines from line number 45500973 to 45500976.


text:
Getting VCF file that contain common SNPs from 6 VCF file using isec
I have 6 VCF files, where I would like to obtain the SNPs that are common (by position) in all the 6 files. I have tried this command
bcftools isec  -p trial2 -n=6 -c all A2_SNP.vcf.gz A_SNP.vcf.gz B_SNP.vcf.gz C_SNP.vcf.gzD_SNP.vcf.gz E_SNP.vcf.gz

But this produces 6 files, and through the README file, I think each of the VCF contains the unique SNPs for each of the 6 files, and I could only get sites.txt file that ocntain the positions that are common. Are there any way to get a VCF file that have the common SNPs ?
Thanks

You've used the right tool and gotten the right results. What you need to understand is that VCF contains variant loci as well as the genotypes for >=1 sample(s) at those loci. The intersect operation will extract common loci, but it does not know which genotypes you're interested in, which is why it generates one file per input, giving you the genotypes from each VCF for the common loci.
You could use the sites file to do a bcftools merge operation, limiting to just the common sites using the -R <regions-file> option so you'd get a VCF with all the samples across all 6 input VCFs but only for the variant loci.

Thanks a lot. For my understanding: Do you means that the output 6 files from isec tool I used is a VCF for each sample and contains the genotype (and other info) for only that common sites (that are also in site.txt) ? If so, do you assume that if I merge this using -R, then that merged VCF file contains all 6 samples. How can I get the genotype or alleles for the 6 individuals ?

Yes, each VCF output should contain information from its corresponding VCF input. You could either use sites file to merge all 6 original inputs or simply merge all 6 isec outputs get to the merged common-loci VCF.

So, merging the 6 VCF file will justput them side by side on  sample basis ? Could you advise on a code for that ?

Try the command, look at the output and come back if you have further questions.


text:
Calculation of Site Frequency Spectrum (SFS) using STR data
Hi everyone I had a question on if it is possible to calculate (Site Frequency Spectrum)SFS using STR data. Please mention if you know any programm that can be used for that!
I have tried Arlequin but it doesn't offer this option for STR. Personally I haven't find anything that can do this type of calculation with STR!
Thank you!

@M__ "I have tried Arlequin but it doesn't offer this option for STR"

There are two things. First STS - which is a complex algorithm, secondly is the data STR. If Arlequin has the maximum likelihood estimate for STS (which the OP  implying it does), it's easy. This is because Arlequin will process any data format available, STR is simply a data format, not an algorithm. It's about looking at what format would be most suitable and converting the data accordingly to a format Arlequin uses. The Arlequin data types are not rigid like phylogenetics. That is a fairly easy question. If Arlequin doesn't have STS (the algorithm), finding an STS algorithms is not trivial.

The parsing operation might be a little complicated - but it's not hard and someone has probably written a parsing script, not necessarily STR -> Arlequin, but STR to an intermediate format that will pass into the official Arlequin data parser.


text:
How Can You Reformat an OTU Table from qiime2 pathway to Include Sequence IDs?
I am trying to run Tax4Fun in MicrobiomeAnalyst, but am running into trouble with my table formatting. From the qiime2 pathway, I have an OTU, metadata and taxonomic table formatted like so:

In order to use Tax4Fun, I understand that I will need to reformat my OTU table to include sequence IDs so it looks like this example:

However, I cannot figure out how to merge my sequence files from the qiime2 pathway with the OTU table so it looks like that. Any suggestions would be much appreciated!

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.


text:
How to handle technical replicates in cellranger count?
I'm working with cellranger count to annotate some publicly available fastq files (SRX7888073). For each sample, 4 runs are available, which I assume represent technical replicates.
How do you recommend to handle them in cellranger? Are they 4 different samples (i.e. 4 separate runs), or 1 sample with 4 lanes?
Thank you !

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.

Welcome @David, nice you joined us. It could be helpful if you could point to the information on the site, ie. "Male 3; Schistosoma mansoni; RNA-Seq". Thats kinda key, clicking offsite links to answer questions is sort of discouraged. The actual experimental details will be explained on the site ...

It depends what you mean by "technical replicates".  If it is the exact same library prep spread out over two lanes, or two flowcells, they must be combined at the fastq stage.  Cellranger will naturally do this with fastqs that keep their original Illumina-assigned names, but just have the files split by lane.
If one file has a read from a gene with a particular cell barcode and UMI, and another read in a separate file but the same library prep has a read from the same gene, with the same barcode and same UMI, those two must not be both counted. So if you combine the counts, that's wrong.
A NextSeq naturally has 4 lanes, so that's almost certainly what you are looking at.  The read names would confirm that, unless they've been altered.  Cat them all together, or alter their names so that they match the original Illumina specs with lanes 1-4.  You might have to do that anyway, cellranger used to be very picky about that, and might still be.

Agreed on all of that. If you alter names though then furst make a new folder, symlink files into that and alter the symlink names. Don't  alter names of actual files, that's error-prone.

If they are truly technical replicates in the most specific sense (e.g. from the same extraction and library prep, just split up for sequencing), then the counts for each identified cell can be combined (i.e. added together). One way to do this is to concatenate the files prior to running any data processing scripts.
If there is anything different about the runs, then they should be treated as non-technical replicates and properly integrated:
https://satijalab.org/seurat/articles/integration_rpca.html
We had a situation with Rhapsody single cell data where a library was split into three on three different cell barcode plates, and was showing a small amount of library preparation bias that influenced clustering. Carrying out an integration of the libraries helped to reduce this bias.

Catting the fastqs will not be the same as combining the counts.  If the catted fastqs are processed by cellranger, deduplication will be handled properly.  If the fastqs are processed separately, you will overcount, because duplicates between runs will not be accounted for.

After some digging, it seems like each of the 4 runs has I1, R1 and R2 files - they're replicates.
You can either

run them separately and then perform cellranger aggr; or
link to all 12 files using soft-links that follow a proper naming convention (the files need to have the same prefix, lane numbers need to be different per replicate), ensure all softlinks are in the same directory, then use the directory with cellranger count.

I'd recommend the former approach as it will probably save a little compute time, is parallelize-able and will be easier to debug/drill down into.

using cell aggr is totally wrong  If you had 4 different samples, it would be a coincidence if some of them happened to have the same barcode, and the right thing to do would be to treat those as two separate cells.  With this, if the same cell barcode turns up in each lane, it's actually the same cell, and they should be combined.


text:
Nextflow - Process has already been used
I am trying to write a nextflow script below:
/*
 * pipeline input parameters
 */
params.reads = "/path/to/fq/files/"
params.guide_library = "/path/to/guidelibrary/"
params.outdir = "results"

log.info """\
    DUAL GUIDE CRISPR QC - N F   P I P E L I N E
    ===================================
    guide library: ${params.guide_library}
reads        : ${params.reads}
    outdir       : ${params.outdir}
    """
    .stripIndent(true)

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

process READCOUNT {
    input:
    path guide_library
    tuple val(sample_name) path(reads)

    output:
    path "$sample_name"

    script:
    """
    exactcounts.py -inputlibrary $guide_library -read1 ${reads[0]} -read2 ${reads[1]} -prefix $sample_name
    """
}

workflow {
    Channel
        .fromFilePairs(params.reads, checkIfExists: true)
        .set { read_pairs_ch }
   count_ch = READCOUNT(params.guide_library,read_pairs_ch)

}

But am getting the following error:
Process 'READCOUNT' has been already used -- If you need to reuse the same component, include it with a different name or include it in a different workflow context

I believe adding a comma after val(sample_name) on line tuple val(sample_name) path(reads) will solve the issue.
EDIT: Also documented here: https://github.com/nextflow-io/nextflow/discussions/2714

Good eye! I had the same and I was sure something was odd but couldn't "see" the missing comma. The resulting error message is really misleading in this case...


text:
Run gmt format download from GSEA
I built the matrix as the gmt example in https://docs.gsea-msigdb.org/#GSEA/Data_Formats/ but It still has errors when I try to run it.
When I try the gene set from https://www.gsea-msigdb.org/gsea/msigdb/human/geneset/HALLMARK_ADIPOGENESIS.html
in .gmt format and load it with GSEA appears this error:
Thank you so much
Nerea

---- Full Error Message ----
col:34 > matrix's fColCnt:34

---- Stack Trace ----
# of exceptions: 1
------col:34 > matrix's fColCnt:34------
java.lang.ArrayIndexOutOfBoundsException: col:34 > matrix's fColCnt:34
at org.gsea_msigdb.gsea/edu.mit.broad.genome.math.Matrix.getColumnV(Matrix.java:261)
at org.gsea_msigdb.gsea/edu.mit.broad.genome.objects.DefaultDataset.getColumn(DefaultDataset.java:289)
at org.gsea_msigdb.gsea/edu.mit.broad.genome.objects.TemplateFactory.extract(TemplateFactory.java:85)
at org.gsea_msigdb.gsea/edu.mit.broad.genome.alg.DatasetGenerators.extract(DatasetGenerators.java:292)
at org.gsea_msigdb.gsea/edu.mit.broad.genome.alg.DatasetGenerators.extract(DatasetGenerators.java:288)
at org.gsea_msigdb.gsea/xtools.gsea.AbstractGsea2Tool.execute_one(AbstractGsea2Tool.java:79)
at org.gsea_msigdb.gsea/xtools.gsea.AbstractGsea2Tool.execute_one_with_reporting(AbstractGsea2Tool.java:103)
at org.gsea_msigdb.gsea/xtools.gsea.Gsea.execute(Gsea.java:166)
at org.gsea_msigdb.gsea/edu.mit.broad.xbench.tui.TaskManager$ToolRunnable.run(TaskManager.java:391)
at java.base/java.lang.Thread.run(Unknown Source)

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.

Please show the first few lines of your input file; from the error, this looks like it might be a formatting issue, and without more information about the format of your file we can't help out.

Hello It was a formatting issue. Sorry, I solved the problem. Thank you

Great! glad it's resolved. Could you kindly add your resolution as an answer please? Firstly it helps the site because unresolved question continue to circulate, secondly it helps your reputation on the site. This is useful should you wish to ask another question here in the future.


text:
Good practices for doing sanity checks when parsing multiple PDB files
This is my first time posting to stack exchange, as well as to using data-driven approaches, so please correct me if my questions is not clear enough.
I am trying to parse 200+ .pdb files to extract residues within 15 angstroms of a certain ligand. I have managed to extract the resnames in a text file.
In order to do a sanity check, I first ran my code on the first 5 pdb files of this list and compared the outputs to expected output. However, I was wondering if there are any other steps I can take to make my check more robust.
Are there any 'gold standard'/good  practices related to this?

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.

Welcome to the site. Could you kindly provide further details, e.g. are you using pymol, biopython PDB and what code have you written to date? Any reproducible example is always welcome.

Add a counter to be sure you're parsing the exact amount of pdbs you expect to parse ? Check for the ligand before extracting the interacting residue ? Check for duplicate pdbs ?

Thank you @M__. I am using VMD to parse through the PDB files. I can share my algorithm as follows:

Load the PDB into vmd. \

Identify the number of occurances of ligand, say Na+. \

For each of the occurance of Na+, identify the atoms within 5A of distance from it. \

Extract the resid and chain names for each of those atoms and store it as a pair in a variable pair_list. \

Filter out the repeated pairs to make sure I am counting every residue only once bold(quick Q: do you think the combo of resid and chain are sufficient to avoid repeated counting?) bold

Extract the resnames of those unique residues and save it in a text file.

Repeat this process for all the PDB in current directory.

I have added the TCL script doing this for NH4+ at this github link: https://github.com/priyanshurg/TCLScripts/blob/main/reslist4-Within5.tcl

Thank you @pippo1980ONSTRIKE!
Could you recommend what property in the PDB file I can use to classify PDBs as similar? For example there are many versions of the KcsA protein in RCSB.
Your point about duplicate pdbs made me realise that my database may be biased by proteins that are studied more often compared to others. Do you have any thoughts on how you'd go about mitigating this bias?

dont know see  https://actanaturae.ru/2075-8251/article/view/10699/pdf materials % methods, you could drop structures too similar or sequences too similar I believe its on you , try both and look at the differences  or you could just concentrate on the active site alone using similar approach