text:
Ensembl Variant Effect Predictor (VEP) issue during execution
I got this error during execution, pls guide me.
vep --cache -i xxxx.vcf -o xxxxx_new.vcf

Possible precedence issue with control flow operator at /software/xxxx_pkgs/envs/ensemble-vep/lib/site_perl/5.26.2/Bio/DB/IndexedBase.pm line 805.
DBI connect('host=ensembldb.ensembl.org;port=3306','anonymous',...) failed: Can't connect to MySQL server on 'ensembldb.ensembl.org' (110) at /software/xxxx_pkgs/envs/ensemble-vep/share/ensembl-vep-88.9-1/Bio/EnsEMBL/Registry.pm line 1771.

-------------------- EXCEPTION --------------------
MSG: Cannot connect to the Ensembl MySQL server at ensembldb.ensembl.org:3306; check your settings & DBI error message: Can't connect to MySQL server on 'ensembldb.ensembl.org' (110)
STACK Bio::EnsEMBL::Registry::load_registry_from_db /software/xxxx_pkgs/envs/ensemble-vep/share/ensembl-vep-88.9-1/Bio/EnsEMBL/Registry.pm:1771
STACK Bio::EnsEMBL::VEP::BaseVEP::registry /software/xxxx_pkgs/envs/ensemble-vep/share/ensembl-vep-88.9-1/modules/Bio/EnsEMBL/VEP/BaseVEP.pm:163
STACK Bio::EnsEMBL::VEP::BaseRunner::setup_db_connection /software/xxxx_pkgs/envs/ensemble-vep/share/ensembl-vep-88.9-1/modules/Bio/EnsEMBL/VEP/BaseRunner.pm:80
STACK Bio::EnsEMBL::VEP::Runner::init /software/xxxx_pkgs/envs/ensemble-vep/share/ensembl-vep-88.9-1/modules/Bio/EnsEMBL/VEP/Runner.pm:72
STACK Bio::EnsEMBL::VEP::Runner::run /software/xxxx_pkgs/envs/ensemble-vep/share/ensembl-vep-88.9-1/modules/Bio/EnsEMBL/VEP/Runner.pm:111
STACK toplevel /software//envs/ensemble-vep/bin/vep:209
Ensembl API version = 88
```

You might prefer to run VEP in offline mode:

Using --cache (without --offline) uses the local cache on disk to
fetch most annotations, but allows database connections for some
features (see cache limitations)

If you need HGVS names in your output, you'll also need to supply a FASTA file using --fasta:

--offline
Enable offline mode. No database connections will be made, and a cache file or GFF/GTF file is required for annotation. Add
--refseq to use the refseq cache (if installed). Not used by default.

--fasta [file|dir]
Specify a FASTA file or a directory containing FASTA files to use to look up reference sequence. The first time you run VEP with this parameter an index will be built which can take a few minutes. This is required if fetching HGVS annotations (--hgvs) or checking reference sequences (--check_ref) in offline mode (--offline), and optional with some performance increase in cache mode (--cache). See documentation for more details. Not used by default.


text:
SyntaxError after end of Snakefile code?
I am trying to run Fastqc on multiple files using Snakemake but am receiving an ambiguous SyntaxError past the last line of code in the Snakefile.
Snakefile:
import pandas as pd

SAMPLE_TABLE=pd.read_csv('./SampleTable.txt',sep='\t')
SAMPLES=SAMPLE_TABLE["sample_name"]
SAMPLES=list(SAMPLES)
DIRECTORY=list(SAMPLE_TABLE["FolderName"])

rule run_fastqc:
    input:
        expand("/PathToInputsFolder/{dir}/{sample}_{{Rnumber}}_001.fastq.gz", zip, dir=DIRECTORY, sample=SAMPLES)
    output:
        expand("{sample}_{{Rnumber}}_001.fastqc.html", sample=SAMPLES)
        expand("{sample}_{{Rnumber}}_001.fastqc.zip", sample=SAMPLES)
    shell:
        "fastqc {input} --outdir=./qc/"

Code:
snakemake --cores 1 {.+}.fastqc.html {.+}.fastqc.zip

Output:
SyntaxError in line 16 of /PathToSnakefile/Snakefile:
invalid syntax (Snakefile, line 16)

How should I figure out where the error is? Thank you for your time!

How to figure out what error messages are trying to tell you is one of the great mysteries… :-)
What version of Snakemake are you using? I get slightly different messages. With an ancient version (3.13.3) on one machine:
SyntaxError in line 14 of /path/to/Snakefile:
invalid syntax

With a more up-to-date version (6.9.1, running on Python 3.9):
SyntaxError in line 16 of /path/to/Snakefile:
invalid syntax. Perhaps you forgot a comma? (Snakefile, line 16)

Both of these are useful hints. So one piece of advice would be to update your Snakemake and/or Python if it is old.
Another useful technique when faced with a mysterious error like this is to temporarily remove parts of your Snakefile until the error goes away. The complicated part of your script towards the end of the file is the output: stanza, so you might try removing one of the expand(…) lines and seeing what happens. If the error goes away, then you have localised the problem and can put the line back and carefully check the syntax in that particular part of the Snakefile.

You're missing a comma between the two output lines.  (John Marshall's more detailed syntax error message was spot-on.)
With a simplified version of your Snakefile:
SAMPLES=["s1", "s2", "s3"]
rule run_fastqc:
    output:
        expand("{sample}_{{Rnumber}}_001.fastqc.html", sample=SAMPLES)
        expand("{sample}_{{Rnumber}}_001.fastqc.zip", sample=SAMPLES)
    shell:
        "echo {output}"

You get the syntax error:
$ snakemake -nrp s1_R1_001.fastqc.html
SyntaxError in line 6 of /data/home/jesse/Snakefile:
invalid syntax (Snakefile, line 6)

With the comma:
SAMPLES=["s1", "s2", "s3"]
rule run_fastqc:
    output:
        expand("{sample}_{{Rnumber}}_001.fastqc.html", sample=SAMPLES), # <-- here
        expand("{sample}_{{Rnumber}}_001.fastqc.zip", sample=SAMPLES)
    shell:
        "echo {output}"

Now it's OK:
$ snakemake -nrp s1_R1_001.fastqc.html
Building DAG of jobs...
Job stats:
job           count    min threads    max threads
----------  -------  -------------  -------------
run_fastqc        1              1              1
total             1              1              1

[Fri Nov 12 08:54:54 2021]
rule run_fastqc:
    output: s1_R1_001.fastqc.html, s2_R1_001.fastqc.html, s3_R1_001.fastqc.html, s1_R1_001.fastqc.zip, s2_R1_001.fastqc.zip, s3_R1_001.fastqc.zip
    jobid: 0
    reason: Missing output files: s1_R1_001.fastqc.html
    wildcards: Rnumber=R1
    resources: tmpdir=/tmp

echo s1_R1_001.fastqc.html s2_R1_001.fastqc.html s3_R1_001.fastqc.html s1_R1_001.fastqc.zip s2_R1_001.fastqc.zip s3_R1_001.fastqc.zip
Job stats:
job           count    min threads    max threads
----------  -------  -------------  -------------
run_fastqc        1              1              1
total             1              1              1

This was a dry-run (flag -n). The order of jobs does not reflect the order of execution.

Each expand() makes a list, so at first Snakemake sees two lists back-to-back and doesn't know what to do with them.  If you use a comma and give output names you could handle the two lists separately in your rule, but if you don't need to, you could also just us a single expand call, like:
output:
    expand("{sample}_{{Rnumber}}_001.fastqc.{ext}", sample=SAMPLES, ext=["html", "zip"])

Other suggestions for troubleshooting:

If the above simplified Snakefile works but the real one still doesn't, edit them incrementally to be more and more similar and you'll hone in on whatever problems are left.
I find trimming away complexity until the problem reveals itself is often easiest, so I'd keep simplifying the real case until that happens.
expand is just a Python function that Snakemake provides you, so you can call it yourself outside of a rule definition and check what it returns.
You can even comment out all the Snakemake-specific syntax and leave just Python code when debugging, to sanity-check the rest of the code (e.g. your inputs and outputs from expand).

You can give someone a fish [a comma] today… or you can teach them to fish [to decrypt their own error messages]…

With the suggested comma, I get the message: "WorkflowError:
Duplicate output file pattern in rule run_fastqc. First two duplicate for entries 0 and 1
File "/PathToSnakefile/Snakefile", line 22, in ". I'm not sure how the output files are duplicate if one is .html and the other is .zip. I also added 'READ=list(SAMPLE_TABLE["read"])' and changed '{{Rnumber}}' to 'R{val}, val=READ' so the wildcard is defined, and it still leads to the duplicate output error :/

I added a few bullets on the end of my answer that might help you.  To John Marshall's point, Snakemake is correctly telling you the problem here, but you're assuming the output files are what you expect them to be.  If you investigate the expand behavior in your new version you'll see the issue.

source: Stack Exchange

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

sort: votes


text:
Gene expression changes in one condition
I want to look at an RNA-seq data-set on 1 condition to look for expressed genes. As there is only one condition, I can’t perform differential analysis.
The only thing I can think to do is set a threshold (i.e. all replicates must have > 10 gene counts) then say these are expressed. But I don’t think this is a very good way of doing it.
Do you know of any better ways that we can do this?

If you've got enough replicates (ideally at least 6, but can be done with 3), you should be able to do a Variance-stabilised Transformation using DESeq2:
vsd <- vst(dds, blind=FALSE)
rld <- rlog(dds, blind=FALSE)
head(assay(vsd), 3)

More details on this can be found in the DESeq2 Vignette:
https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#extracting-transformed-values

Answer from @james-hawley, converted from comment:
In any high throughput sequencing experiment, there will be lowly expressed RNA molecules that are not captured by the library prep + sequencing. So defining what "is expressed" is somewhat arbitrary from RNA-seq data, regardless of how many replicates you have and what method you use. That said, you can control some of the noise in your RNA-seq data using a variance-stabilized transformation, like @gringer suggests, and look at the distribution of genes/transcripts. Some threshold like $rlog \ge x$ in at least $n$ samples is usually a good starting point (for some arbitrary $x$ and $n$).

I don't think there is a very good way to do what you want to do.  I'd think some kind of RPM would be better than strictly raw counts; you have to account for overall library size.  Or TPM might be better.

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


text:
Finding bioperl on cpanmetadb failed
I failed to install bioperl on debian in the following way:
$ sudo cpanm bioperl
! Finding bioperl on cpanmetadb failed.
! Finding bioperl () on mirror http://www.cpan.org failed.
! Couldn't find module or a distribution bioperl

Thank you in advance,

conda install -c bioconda perl-bioperl  instead of using CPAN


text:
xtract from eutilities not exctracting the information I need
I'm running the following command which uses the NCBI e-utilities:
esearch -db bioproject -query "PRJNA198476 [PRJA]" | esummary
This returns:
<?xml version="1.0" encoding="UTF-8" ?>
<!DOCTYPE DocumentSummarySet PUBLIC "-//NLM//DTD esummary bioproject 20140903//EN" "https://eutils.ncbi.nlm.nih.gov/eutils/dtd/20140903/esummary_bioproject.dtd">

<DocumentSummarySet status="OK">
<DbBuild>Build211111-1430.1</DbBuild>

<DocumentSummary><Id>198476</Id>
    <TaxId>6239</TaxId>
    <Project_Id>198476</Project_Id>
    <Project_Acc>PRJNA198476</Project_Acc>
    <Project_Type>Primary submission</Project_Type>
    <Project_Data_Type>Transcriptome or Gene expression</Project_Data_Type>
    <Sort_By_ProjectType>476331</Sort_By_ProjectType>
    <Sort_By_DataType>472932</Sort_By_DataType>
    <Sort_By_Organism>51608</Sort_By_Organism>
    <Project_Subtype></Project_Subtype>
    <Project_Target_Scope>Multiisolate</Project_Target_Scope>
    <Project_Target_Material>Transcriptome</Project_Target_Material>
    <Project_Target_Capture>Whole</Project_Target_Capture>
    <Project_MethodType>Sequencing</Project_MethodType>
    <Project_Method></Project_Method>
    <Project_Objectives_List>
        <Project_Objectives_Struct>
            <Project_ObjectivesType>Expression</Project_ObjectivesType>
            <Project_Objectives></Project_Objectives>
        </Project_Objectives_Struct>
    </Project_Objectives_List>
    <Registration_Date>2013/04/22 00:00</Registration_Date>
    <Project_Name>Caenorhabditis elegans</Project_Name>
    <Project_Title>Deep sequencing of endogenous mRNA from Caenorhabditis elegans in the presence and absence of arsenite</Project_Title>
    <Project_Description>Background: Arsenite is one of the most toxic chemical substances known and is assumed to exert detrimental effects on viability even at lowest concentrations. By contrast and unlike higher concentrations, we here find that exposure to low-dose arsenite promotes growth of cultured mammalian cells. In the nematode C. elegans, low-dose arsenite promotes resistance against thermal and chemical stressors, and extends lifespan of this metazoan, whereas higher concentrations reduce longevity. While arsenite causes a transient increase in reactive oxygen species (ROS) levels in C. elegans, co-exposure to ROS scavengers prevents the lifespan-extending capabilities of arsenite, indicating that transiently increased ROS levels act as transducers of arsenite effects on lifespan, a process known as mitohormesis. The RNA-seq data comprises 2 biological replicates for worms exposed to 100nM Arsenite 48h after L4 and 2 biological replicates of the same age as controls Jena Centre for Systems Biology of Ageing - JenAge (www.jenage.de) Overall design: 4 samples: 2 mRNA profiles of C.elegans 48h after L4 exposed to Arsenite; 2 mRNA profiles of C.elegans 48h after L4 as controls (H20). The N2 wild type (var. Bristol) strain was used.</Project_Description>
    <Keyword></Keyword>
    <Relevance_Agricultural></Relevance_Agricultural>
    <Relevance_Medical></Relevance_Medical>
    <Relevance_Industrial></Relevance_Industrial>
    <Relevance_Environmental></Relevance_Environmental>
    <Relevance_Evolution></Relevance_Evolution>
    <Relevance_Model>Yes</Relevance_Model>
    <Relevance_Other></Relevance_Other>
    <Organism_Name>Caenorhabditis elegans</Organism_Name>
    <Organism_Strain></Organism_Strain>
    <Organism_Label></Organism_Label>
    <Sequencing_Status>SRA/Trace</Sequencing_Status>
    <Submitter_Organization>Leibniz Institute for Age Research - Fritz Lipmann Institute</Submitter_Organization>
    <Submitter_Organization_List>
        <string>Leibniz Institute for Age Research - Fritz Lipmann Institute</string>
    </Submitter_Organization_List>
    <Supergroup>Eukaryotes</Supergroup>
</DocumentSummary>

</DocumentSummarySet>

I'm trying to extract the Submitter_Organization using:
esearch -db bioproject -query "PRJNA198476 [PRJA]" | esummary | xtract -pattern DocummentSummary -element Submitter_Organization
However, this returns nothing, and I can't figure out why

You have a typo in DocumentSummary.

@vkkodali thank you that was it


text:
Does anyone know what tools are used to construct this beautiful genome map?

Does anyone know what tools are used to construct this beautiful genome map?
The research paper didn't mention how to construct this map.
Link to paper: https://www.sciencedirect.com/science/article/pii/S2052297515000773?via%3Dihub

Circos: http://circos.ca/ (or a copy of it written for another language, depending on what the paper's authors used — probably best to check with them).

hey thanks! any tutorial on how to make those maps in circos? i found out that circos kinda hard to use for beginners

To expand on this, Circos is an older Perl executable script, a mighty complex one with loads of functionality, but awkward input that requires a few trials. For scripting in Python/Julia/JS the Plotly module has a tool, albeit it slightly trickier to overlay stuff one than for a scatter, while in R there are a few options. the RCircos is the R package that is more faithful to the original Perl functionality I believe.

The quick start guides are the best way to get going with the original Circos: circos.ca/documentation/tutorials/quick_guide It looks complicated, but it is no more difficult to use than any of its imitations, which all require as much work in the end, as the same data are required to make the desired figure.

I tried to used Circos for my MSc Thesis but didn't have enough time to actually learn how to use so I managed to find a less fancy but also cool version of it called ClicO.FS which is an interactive web-based service of Circos (as the title of the article below says).

Possibly BRIG (blast ring image generator) was used
http://brig.sourceforge.net/
I also catalog other tools here if interested https://cmdcolin.github.io/awesome-genome-visualization/?latest=true&selected=%23BRIG

Circos appears to have been published at least two years before BRIG, but no citation. Unfortunate.


text:
Download COX1 (COI) gene via biopython using accessions for entire mitochondrial genomes
I have a list of accessions for the the entire mitochondrial genomes for big cats. I need to download the COX1 genes for each of these accessions.
Here is one accession and here is a link to its COX1 gene, which I found manually on that page.
I have downloaded all the accessions as an XML genbank file:
#!/usr/bin/env python3

from Bio import Entrez
from Bio import GenBank
from Bio import SeqIO
import os

if os.path.exists('cats.xml') ==False:
    with open('accessions.txt', 'r') as fin:
        #for line
        accs = [line.rstrip() for line in fin ]

    print(accs)

    Entrez.email =  "first.last@example.se"
    rec = Entrez.efetch(db='nuccore', id=accs, rettype="gb", retmode="xml")
    with open('cats.xml', 'w') as out:
        print(rec.read().decode("utf-8"), file=out)

Then I tried to parse this file:
with open('cats.xml', 'rb') as xin:
    for record  in Entrez.parse(xin):
        #print(dir(record))
        #print(record.keys())
        pass
print(record['GBSeq_feature-table'][0])

But it's a mess of nested lists and dictionaries. Help!

Please [edit] your question and show us the output you expect from your example input. I get what looks like a perfectly normal and valid JSON output. What exactly do you need? COX1 gene accessions? DNA sequence? cDNA sequence? Protein translation?


text:
Error while parsing gene bank file using Biopython
I am trying to extract the protein sequence of specific genes from gene bank like format file obtained from antismash part of which looks like this
LOCUS       scaffold_10            47160 bp    DNA     linear   UNK 01-JAN-1980
DEFINITION  scaffold_10.
ACCESSION   scaffold_10
VERSION     scaffold_10
KEYWORDS    .
SOURCE
  ORGANISM
            .
COMMENT     ##antiSMASH-Data-START##
            Version      :: 6.0.1-25e43a7(changed)
            Run date     :: 2021-11-05 17:39:45
            NOTE: This is a single cluster extracted from a larger record!
            Orig. start  :: 2049841
            Orig. end    :: 2097001
            ##antiSMASH-Data-END##
FEATURES             Location/Qualifiers
     protocluster    1..47160
                     /aStool="rule-based-clusters"
                     /contig_edge="False"
                     /core_location="join{[2069841:2070128](+),
                     [2070175:2070449](+), [2070507:2071973](+),
                     [2072023:2074178](+), [2074232:2074784](+),
                     [2074838:2077001](+)}"
                     /cutoff="20000"
                     /detection_rule="cds(PKS_AT and (PKS_KS or ene_KS or mod_KS
                     or hyb_KS or itr_KS or tra_KS))"
                     /neighbourhood="20000"
                     /product="T1PKS"
                     /protocluster_number="1"
                     /tool="antismash"
     proto_core      join(20001..20287,20335..20608,20667..22132,22183..24337,
                     24392..24943,24998..27160)
                     /aStool="rule-based-clusters"
                     /tool="antismash"
                     /cutoff="20000"
                     /detection_rule="cds(PKS_AT and (PKS_KS or ene_KS or mod_KS
                     or hyb_KS or itr_KS or tra_KS))"
                     /neighbourhood="20000"
                     /product="T1PKS"
                     /protocluster_number="1"
     cand_cluster    1..47160
                     /SMILES="CC(=O)C(=O)O"
                     /candidate_cluster_number="1"
                     /contig_edge="False"
                     /detection_rules="cds(PKS_AT and (PKS_KS or ene_KS or
                     mod_KS or hyb_KS or itr_KS or tra_KS))"
                     /kind="single"
                     /product="T1PKS"
                     /protoclusters="1"
                     /tool="antismash"
     region          1..47160
                     /candidate_cluster_numbers="1"
                     /contig_edge="False"
                     /product="T1PKS"
                     /region_number="1"
                     /rules="cds(PKS_AT and (PKS_KS or ene_KS or mod_KS or
                     hyb_KS or itr_KS or tra_KS))"
                     /tool="antismash"
     CDS             1899..2720
                     /ID="CDS_13110"
                     /NRPS_PKS="Domain: PKS_KR (26-149). E-value: 7.6e-08.
                     Score: 24.2. Matches aSDomain:
                     nrpspksdomains_jgi.p_FusspF11_1_766668_PKS_KR.1"
                     /NRPS_PKS="type: other"
                     /gene="jgi.p_FusspF11_1_766668"
                     /gene_functions="biosynthetic-additional (smcogs)
                     SMCOG1001:short-chain dehydrogenase/reductase SDR (Score:
                     161.6; E-value: 3.9e-49)"
                     /gene_kind="biosynthetic-additional"
                     /phase="0"
                     /source="prediction"
                     /transl_table=1
                     /translation="MTPPPKSALESDPAVGAQRRFSVTGNAVVTGGAGVLGLHACDALL
                     EHGLEGLMILDVNPAQSQGQITSLQNKFPRAKIMALKVDVTDENVVNAAMEETARVLGS
                     IDTLICFVGVVGCVETLDMPVPQWRKILDINTTGSFICAQAAARQMVKRGRGGSIVFVA
                     SISAHRVNYPQPQAAYNVSKSALLMLKSCLAAEWARYGIRTNSISPGYMDTILNEGDGI
                     AEHRKIWAEHNPSGRMGAPSELTGTVVLLASSAGSYINGADIVVDGGGIVL"
     aSDomain        1977..2345
                     /aSDomain="PKS_KR"
                     /aSTool="nrps_pks_domains"
                     /database="nrpspksdomains.hmm"
                     /detection="hmmscan"
                     /domain_id="nrpspksdomains_jgi.p_FusspF11_1_766668_PKS_KR.1
                     "
                     /evalue="7.60E-08"
                     /label="jgi.p_FusspF11_1_766668_PKS_KR.1"
                     /locus_tag="jgi.p_FusspF11_1_766668"
                     /protein_end="149"
                     /protein_start="26"
                     /score="24.2"
                     /specificity="KR activity: inactive"
                     /specificity="KR stereochemistry: C1"
                     /tool="antismash"
                     /translation="AVVTGGAGVLGLHACDALLEHGLEGLMILDVNPAQSQGQITSLQN
                     KFPRAKIMALKVDVTDENVVNAAMEETARVLGSIDTLICFVGVVGCVETLDMPVPQWRK
                     ILDINTTGSFICAQAAARQ"
     aSModule        1977..2345
                     /domains="nrpspksdomains_jgi.p_FusspF11_1_766668_PKS_KR.1"
                     /incomplete
                     /locus_tags="jgi.p_FusspF11_1_766668"
                     /tool="antismash"
                     /type="pks"
     CDS_motif       1980..2042
                     /aSTool="nrps_pks_domains"
                     /database="abmotifs"
                     /detection="hmmscan"
                     /domain_id="nrpspksmotif_jgi.p_FusspF11_1_766668_0001"
                     /evalue="1.90E-05"
                     /label="PKSI-KR_m1"
                     /locus_tag="jgi.p_FusspF11_1_766668"
                     /protein_end="48"
                     /protein_start="27"
                     /score="17.3"
                     /tool="antismash"
                     /translation="VVTGGAGVLGLHACDALLEHG"
     CDS             join(3252..3263,3321..3443,3493..4428)
                     /ID="CDS_13111"
                     /NRPS_PKS="Domain: PKS_ER (28-277). E-value: 6.4e-13.
                     Score: 40.2. Matches aSDomain:
                     nrpspksdomains_jgi.p_FusspF11_1_794272_PKS_ER.1"
                     /NRPS_PKS="type: other"
                     /gene="jgi.p_FusspF11_1_794272"
                     /gene_functions="biosynthetic-additional (smcogs)
                     SMCOG1040:alcohol dehydrogenase (Score: 275.7; E-value:
                     7.6e-84)"
                     /gene_kind="biosynthetic-additional"
                     /phase="0"
                     /source="prediction"
                     /transl_table=1
                     /translation="MSQTNLSCVLYGPGKARFENRPVPSLKDAHDVIIRISYVGVCGSD
                     VHFWTDGGFARKVSEDQPLVMGHEASGIVRSIGPDVTLLKPGDRVAIEPGFSCRRCKQC
                     KDGRYNLCPKMKFAADPPLTQGTLSRFFSIPEDFAYKIPDSLSLEEAVLVEPLAVAVHG
                     IRLAGLEVGQRVLVQGSGTIGLLTAAVAKAYGAKQVYITDVNLDKIKFAKKYLECSAFI
                     PDLGSTPEENAARFKTETGLDDGVDAVIECTGVEASTQTGLLALSAGGVLVQVGLGKPV
                     QAIPIHAMSEKEIVLKTSFRYGPGDYEIALELLESGKVSVRPLISSITPFEKATEAWEK
                     TRKGKGIKNLIRGVQD"

Now from this I want to extract the protein sequence of gens having PKS domain (PKS_ER/PKS_AT/PKS_KS ...). So far I have tried the following biopython code but it is giving me error
from Bio import SeqIO
from Bio.SeqFeature import SeqFeature

gbk = "myfile_predicted_SMclusters.gbk"
fa = "myfile_PKS_protein.fa"
input_handle = open(gbk, "r")
output_handle = open(fa, "w")

for record in SeqIO.parse(input_handle, "genbank"):
    features = [feature for feature in record.features if feature.type == "CDS"]
    for feature in features:
        if feature.qualifiers["NRPS_PKS"][1] == "PKS":
        assert len(seq_feature.qualifiers['translation'])==1
        output_handle.write(">%s %s\n%s\n" % (
            record.description,
            feature.location,
            seq_feature.qualifiers['translation'][0]))

output_handle.close()
input_handle.close()

The error is as follows:
if feature.qualifiers["NRPS_PKS"][1] == "PKS":
KeyError: 'NRPS_PKS'

The link to the sample inputfile is below:
sample my_file.gbk
I'm using python 3.8.8
Please help!!!

I see 2 problems.
First, you need to check if the qualifier your are searching for exists in the given CDS feature.
Second, the qualifier holds more then just the string PKS for example Domain: PKS_AT (537-846). E-value: 1.5e-98.
You could modify your code like this:
feature_type = "CDS"
feature_qualifier = "NRPS_PKS"
target = "PKS"

def has_pks_value(feature, qualifier, part):
    return any(part in target for target in feature.qualifiers[qualifier])

for record in SeqIO.parse(input_handle, "genbank"):
    for feature in record.features:
        if feature.type == feature_type and feature_qualifier in feature.qualifiers.keys() and has_pks_value(feature, feature_qualifier, target) and len(feature.qualifiers['translation'])==1:
            output_handle.write(f">{record.description} {feature.location.start} {feature.location.end} {feature.location.strand} \n {SeqFeature(feature.location).extract(record.seq)}")

I've created some variables which for your string. This might prevent you from making some typos.
Also instead of a list comprehension, I just run over all features and checked with a single if statement if the asked qualifier exists and contains the given value.
I also changed the string formatting to f-string. But I think you should create the Fasta Files using the SeqIO module of Biopython.
As result I got 37 sequences

Update - I added a missing method to check for the correct value of the qualifier

Thank you for your suggestions, But I god the following error now.  output_handle.write(f">{record.description} {feature.location.start} {feature.location.end} {feature.location.strand} \n {SeqFeature(feature.location).extract(record.seq)}")                                                                                ^ SyntaxError: invalid syntax
could you please Help, I recently started using python (last 3 days) just for this purpose only trying to make sense of resources available online.

Could you try run it without the f-string part? Just use your original version for the file writing

It has another error output_handle.write(">%s %s\n%s\n" % (record.description, feature.location, seq_feature.qualifiers['translation'][0])) NameError: name 'seq_feature' is not defined

The variable seq_feature doesn't exist anymore. Use feature instead


text:
How to map list items which are in groups in the dictionary values and return keys corrosponding it?
dict1 = {'chr1:1000-2000': ['A', 'B', 'C'], 'chr2:999-2000': ['A', 'X', 'Y'], 'chr3:1200-3000': ['B', 'A', 'C']}

list1 = [('A', 'B'), ('A', 'C'), ('A', 'X'), ('A', 'Y'), ('B', 'C'), ('B', 'X'), ('B', 'Y'), ('C', 'X'), ('C', 'Y'), ('X', 'Y')]

I want to map the items in list1 eg('A','B') both in dict1 values and return keys in wich this group occur

Please [edit] your question and show us the expected output. Should both keys be there or is one enough? Can they be in any order?

edited thank you

I am assuming you only want keys where both tuple items are present. If so, this should work:
no_match = []
for item in list1:
    item1, item2 = item
    wanted_keys = [
        key for key in dict1.keys() if
        item1 in dict1[key] and
        item2 in dict1[key]
    ]
    if len(wanted_keys) > 0:
        print("%s:  %s" % (item, wanted_keys))
    else:
        no_match.append(item)
print("No matches found for items:", no_match)

The output looks like this:
('A', 'B'):  ['chr1:1000-2000', 'chr3:1200-3000']
('A', 'C'):  ['chr1:1000-2000', 'chr3:1200-3000']
('A', 'X'):  ['chr2:999-2000']
('A', 'Y'):  ['chr2:999-2000']
('B', 'C'):  ['chr1:1000-2000', 'chr3:1200-3000']
('X', 'Y'):  ['chr2:999-2000']
No matches found for items: [('B', 'X'), ('B', 'Y'), ('C', 'X'), ('C', 'Y')]

There is probably a better, more elegant way of doing this, I am not a python person, but this is the best I could come up with.


text:
Synteny test and tools to use for mammalian organisms?
I am trying to do a synteny test for MHC genes in human and another mammalian organism
For this reason, I need to download all paralogs of rabbits from humans and then construct position the coordinates of paralog to test the synteny. However, I am unsure how to download the paralogs from Ensembl and how to perform this test?
My main aim is to research synteny between MHC genes of humans, rabbits and mice.
Any suggestions on how to do this or alternative bioinformatic methods on how to perform a synteny test would be highly appreciated!

Are you sure you mean paralogues (genes separated by a duplication event) and not orthologues (genes separated by a speciation event)? Please confirm this. I will separately answer your question assuming you meant orthologues.

This is not an easy ask. Firstly, very few species (I think only human and mouse) have alternative haplotypes annotated for MHC. Most species only have the primary assembly, only one possible haplotype of this locus. This means that the haplotype in the primary assembly is not necessarily the same haplotype as you'll find in another species. Since the MHC haplotypes can contain completely different genes, you will necessarily be able to find orthologues and synteny.
Secondly, we in Ensembl do not calculate homologues of the genes on the alternative loci, so we only have homologues of the genes that just happen to be on the primary assembly.
There is already pre-calculated synteny in Ensembl. Here's the human MHC locus compared to Rabbit. This includes all the homologous genes across this locus. However, this is just one arbitrary MHC vs another arbitrary MHC.


text:
Peptide data visualization
I'm trying to replicate the figure given in this paper using this table which they have given.

So what i understand from the figure is they have calculated the zscore between two condition and plotted the same.
Now How is the length of each bar which is the peptide here is variable is it the zscore or the actual peptide length?
Small subset of the data
dput(head(d2))
structure(list(Peptide = structure(1:6, .Label = c("AAHLPAEFTPAVH",
"AAHLPAEFTPAVHASLD", "AAWGKVG", "AEFTPAVHA", "AEFTPAVHASLDKFLASVSTVLTSKYR",
"AFSDGLAH", "AFSDGLAHLD", "AFSDGLAHLDNLK", "AGEYG", "AHAGEYG",
"AHAGEYGAEALE", "AHAGEYGAEALERMF", "AHAGEYGAEALERMFLSFPTTK",
"AHGKKVLG", "AHLDNLK", "AHLPAEFTPAVH", "AHLPAEFTPAVHASLD", "AHLPAEFTPAVHASLDKF",
"AHVDDMPNA", "AHVDDMPNALS", "ALSDLH", "ALWGKVNVDEVG", "ALWGKVNVDEVGGEA",
"ALWGKVNVDEVGGEALG", "ALWGKVNVDEVGGEALGRLL", "ANALAHKYH", "ASLDKF",
"ASLDKFLA", "ASLDKFLASVS", "ASLDKFLASVSTVLT", "ASLDKFLASVSTVLTS",
"ASLDKFLASVSTVLTSKYR", "ATLS", "AVAHVDD", "AVAHVDDMPN", "AVAHVDDMPNA",
"AVAHVDDMPNAL", "AVAHVDDMPNALS", "AVAHVDDMPNALSA", "AYQKVVAGVANALAHK",
"DALT", "DALTN", "DDMPNA", "DDMPNALS", "DEVG", "DEVGGEA", "DEVGGEALG",
"DEVGGEALGRLL", "DGLA", "DGLAH", "DGLAHLD", "DGLAHLDNLK", "DKLHVD",
"DKLHVDPEN", "DKLHVDPENF", "DKLHVDPENFR", "DKLHVDPENFRL", "DKLHVDPENFRLL",
"DKLHVDPENFRLLG", "DKTNVK", "DLH", "DLS", "DLSTPD", "DLSTPDA",
"DMPNALS", "DPENF", "DPENFR", "DPVNF", "DPVNFK", "EFTPAVH", "EFTPAVHASLD",
"EFTPAVHASLDKFLASVSTVLTSKYR", "EFTPPVQ", "EFTPPVQA", "EFTPPVQAA",
"ESFG", "ESFGDL", "ESFGDLSTPDAVM", "ESFGDLSTPDAVMG", "EVGG",
"EVGGEALG", "FDLS", "FESF", "FESFG", "FESFGDLS", "FESFGDLSTPDAVMGNPKVK",
"FGDL", "FGDLS", "FPHF", "FPHFD", "FPHFDL", "FPHFDLS", "FPT",
"FPTT", "FPTTKT", "FPTTKTY", "FPTTKTYFPHF", "FSDGLA", "FSDGLAH",
"FSDGLAHLD", "FSDGLAHLDNLK", "FTPAVH", "FTPAVHASLD", "FTPAVHASLDKFLASVSTVLT",
"FTPPVQ", "FTPPVQA", "FTPPVQAA", "GAFSDGLA", "GAFSDGLAH", "GAFSDGLAHLDNLK",
"GEALG", "GEYG", "GGEAL", "GKKVL", "GKKVLG", "GKVNVD", "GKVNVDE",
"GKVNVDEVG", "GKVNVDEVGGEA", "GKVNVDEVGGEAL", "GKVNVDEVGGEALG",
"GKVNVDEVGGEALGRL", "GKVNVDEVGGEALGRLL", "GLAHLDNLK", "GNPKV",
"GNPKVK", "GNVL", "GRLL", "GTFAT", "GVANALAHKYH", "HAGEYG", "HAGEYGAEALE",
"HAGEYGAEALERMF", "HFDLS", "HHFG", "HKYH", "HLD", "HLDNLK", "HLPAE",
"HLPAEF", "HLPAEFTPAVH", "HLPAEFTPAVHASLD", "HLTPEEK", "HVDD",
"HVDDMP", "HVDDMPN", "HVDDMPNA", "HVDDMPNAL", "HVDDMPNALS", "HVDDMPNALSA",
"HVDDMPNALSALS", "HVDP", "HVDPEN", "HVDPENF", "HVDPENFR", "HVDPENFRL",
"HVDPENFRLLG", "KAAWGKVG", "KEFTPPVQ", "KEFTPPVQA", "KEFTPPVQAA",
"KEFTPPVQAAYQ", "KLHVDPENFRLLG", "KTYFP", "KTYFPH", "KTYFPHF",
"KTYFPHFD", "KVNVDEVGG", "KVNVDEVGGEA", "KVNVDEVGGEALG", "LDNL",
"LDNLK", "LGAFSDGLA", "LHVDP", "LHVDPENF", "LHVDPENFRLLG", "LKGTFATL",
"LLVVYPWTQRF", "LPAE", "LPAEF", "LPAEFTPAVH", "LPAEFTPAVHASLD",
"LPAEFTPAVHASLDK", "LPAEFTPAVHASLDKF", "LSALS", "LSFPTTK", "LSFPTTKT",
"LSFPTTKTY", "LSFPTTKTYFPHF", "LSFPTTKTYFPHFDL", "LSFPTTKTYFPHFDLS",
"LSFPTTKTYFPHFDLSHG", "LSTPDA", "LTPEEK", "LTPEEKSAVT", "LTPEEKSAVTA",
"LTPEEKSAVTAL", "LTPEEKSAVTALW", "LTPEEKSAVTALWG", "LTPEEKSAVTALWGKVNVD",
"LTPEEKSAVTALWGKVNVDEVG", "LTPEEKSAVTALWGKVNVDEVGGEA", "LTPEEKSAVTALWGKVNVDEVGGEALG",
"LVVYPWTQRF", "LVVYPWTQRFF", "LWGKVNVDEVGGEALG", "NAVAHVDD",
"NAVAHVDDMPNA", "NAVAHVDDMPNAL", "NAVAHVDDMPNALS", "NAVAHVDDMPNALSA",
"NAVAHVDDMPNALSALS", "NAVAHVDDMPNALSALSDL", "NAVAHVDDMPNALSALSDLHAHKLR",
"NVDEVG", "NVDEVGG", "NVDEVGGEA", "NVDEVGGEALG", "PADKTNVK",
"PAEF", "PAEFTPAVH", "PAEFTPAVHASLD", "PAVH", "PDAV", "PDAVM",
"PDAVMG", "PEE", "PEEK", "PEEKSAVT", "PEEKSAVTAL", "PENF", "PENFR",
"PENFRLLG", "PHF", "PHFD", "PHFDL", "PHFDLS", "PKVK", "PNAL",
"PNALS", "PPVQ", "PPVQA", "PPVQAA", "PVNF", "PVNFK", "PVQ", "PVQA",
"PWTQ", "PWTQRF", "PWTQRFF", "QVKGHGKKVADALTN", "RVDPVNF", "SDGLA",
"SDGLAHLDN", "SDGLAHLDNLK", "SDLH", "SFGD", "SFGDL", "SFGDLS",
"SFGDLSTPDAVMG", "SLDKF", "SPADK", "SPADKTNVK", "SVSTVLT", "TKTYFPHF",
"TNVK", "TPAVH", "TPAVHASLD", "TPAVHASLDKFLASVSTVLT", "TPDAVM",
"TPDAVMG", "TPDAVMGNPKVK", "TPEEKSAVT", "TPEEKSAVTAL", "TPEEKSAVTALW",
"TPPVQ", "TSKY", "TVLT", "TYFP", "TYFPH", "TYFPHF", "TYFPHFD",
"TYFPHFDLS", "VAHVDD", "VAHVDDMPNA", "VAHVDDMPNALS", "VAHVDDMPNALSA",
"VDDMP", "VDDMPN", "VDDMPNA", "VDDMPNAL", "VDDMPNALS", "VDEVG",
"VDEVGG", "VDEVGGEA", "VDEVGGEALG", "VDPENF", "VDPENFR", "VDPENFRL",
"VDPENFRLL", "VDPENFRLLG", "VDPVN", "VDPVNF", "VDPVNFK", "VDPVNFKL",
"VDPVNFKLL", "VDPVNFKLLS", "VHL", "VHLT", "VHLTP", "VHLTPE",
"VHLTPEE", "VHLTPEEK", "VHLTPEEKSAVT", "VHLTPEEKSAVTAL", "VHLTPEEKSAVTALWGKVNVD",
"VHLTPEEKSAVTALWGKVNVDE", "VHLTPEEKSAVTALWGKVNVDEVG", "VHLTPEEKSAVTALWGKVNVDEVGGEA",
"VHLTPEEKSAVTALWGKVNVDEVGGEALG", "VHLTPEEKSAVTALWGKVNVDEVGGEALGRLL",
"VKAAWGKVGAHAGEY", "VLSP", "VLSPA", "VLSPAD", "VLSPADK", "VLSPADKT",
"VLSPADKTN", "VLSPADKTNVK", "VLSPADKTNVKAA", "VLSPADKTNVKAAWGKVG",
"VNVDE", "VNVDEVG", "VNVDEVGG", "VNVDEVGGEA", "VNVDEVGGEAL",
"VNVDEVGGEALG", "VNVDEVGGEALGRL", "VVAGVANALA", "VVYPWTQRF",
"VVYPWTQRFF", "VVYPWTQRFFE", "VVYPWTQRFFESFGDLS", "VYPWTQRF",
"VYPWTQRFF", "VYPWTQRFFE", "WGKVG", "WGKVNVD", "WGKVNVDE", "WGKVNVDEVG",
"WGKVNVDEVGGEA", "WGKVNVDEVGGEAL", "WGKVNVDEVGGEALG", "WGKVNVDEVGGEALGRL",
"WGKVNVDEVGGEALGRLL", "YFPHFDLS", "YPWT", "YPWTQ", "YPWTQR",
"YPWTQRF", "YPWTQRFF", "YPWTQRFFE"), class = "factor"), C2 = c(-1.13005927266738,
-1.15451123584162, -0.771975424114207, -0.891503793435265, 0.686928914375861,
-1.15469327326147), C4 = c(0.770515719059141, 0.559148858503814,
1.12965461223673, 1.08129555353272, 0.46033682419062, 0.58089396147208
), C6 = c(0.359543553608241, 0.595362377337809, -0.357679188122518,
-0.18979176009746, -1.14726573856648, 0.573799311789392)), row.names = c(NA,
6L), class = "data.frame")

Any suggestion or help would be really appreciated

The length of the bars and their placement on the x-axis represents the length of the peptide and mapping to the source protein, respectively. The coloring represents the average difference in peptide MS intensity between the control and two experiments.

so i need the length information isnt it which is not given in that paper ..

it can be done using ggplot i mean mapping the length of the peptide?

Maybe I am not following your question, but you get the peptide lengths from the peptides sequences. What you are missing from the table you link to is the source protein sequences, HBb and Hba. The placement of the peptide in the source protein is your x value of the bars, and the length of the peptide is the height of the bar.

"What you are missing from the table you link to is the source protein sequences, HBb and Hba." yes this one . "The placement of the peptide in the source protein is your x value of the bars, and the length of the peptide is the height of the bar." color mapping is the zscore which i do have ..


text:
FindMarkers from Seurat returns p values as 0 for highly significant genes
I have been working on FindMarkers function for identifying significant genes in the cluster. But some Significant genes have very low p values, so they are returned as 0 in the output.Any value less than the .machine$double.xmin limit on R(that is 1e-305). So this output cannot be used for further plots like Volcano plots for visualization.
Any solution to get p values for all the genes?
cd4.markers <- FindMarkers(sc.combined, ident.1 = "CD4+ T cell", min.pct = 0.25)

Part of the cd4.markers output

you can add some constant to those columns

Since other values are proportionally related, if we use a constant those genes forms a cluster.   -Thanks

I prefer to apply a threshold when showing Volcano plots, displaying any points with extreme / impossible p-values (e.g. < 10e-20) with a different symbol at the top of the graph.
However, a better approach is to avoid using p-values as quantitative / rankable results in plots; they're not meant to be used in that way. You can show results as a MA plot instead, plotting log2 fold change vs average expression:

It's a bit trickier to get expression values out of Seurat because they're not currently calculated in the FindMarkers results tables, so you'll need to manually subset the cells and calculate mean expression on a per-marker basis.

It was really helpful solution. I used the AverageExpression(object = pbmc_small) to get mean expression per cell type. then combed it with the Find Markers results using the Gene names then was able to perform MA plot. So those genes are visible under the plot now. - Thanks a lot


text:
What exactly is a reasonable evolutionary criterion?
I got a rejection of publication which the reviewers mainly point out as not following a reasonable evolutionary criterion for choice of species used for phylogeny creation. In our dataset, a list of 100+ eukaryotic species was used to avoid bias. What exactly is a good "reasonable evolutionary criterion"? We are studying an ancient protein (receptor).phylogeny is created based on protein sequences.Can someone provide insight into this? Any suggestions would be appreciated.


text:
How to analyze an IGV alignment
I'm working on a project where I am analyzing the performance of an alignment workflow. My goal is to find regions in the resulting BAM file where there are outstanding discrepancies or anything that indicates my assembly/alignment has "mistakes". I'm working with human genomes from the 1000 genomes project, they are high coverage Illumina paired end reads, each read is 100bp.
My workflow so far:

input paired end FASTQ files (human) into SPAdes to create contigs.fasta file
align paired end FASTQ files to contigs.fasta file with minimap2 to get SAM file
convert SAM to BAM with samtools
sort and index BAM file
load BAM file into IGV and load contigs.fasta as the reference

In the IGV GUI, I have colored alignments by insert size and pair orientation. This has revealed a bunch of insertions, tandem duplications, and inversions.
I am wondering:

Do these colored reads indicate that the assembler made "mistakes" or could they be there just by chance?
Do SNPs indicate anything about how SPAdes/minimap2/samtools performed?
Is there a better way to sort/color the alignments based on what I'm trying to do?

Thank you so much to anyone who has any advice!
EDIT: here is a screenshot of the alignment. For clarification, there are 8 genomes of interest, hence 8 tracks

Background
Integrative genome viewer (IGV) is a well known Broad Institute tool and available here.
The most recent papers are:
James T. Robinson, Helga Thorvaldsdóttir, Aaron M. Wenger, Ahmet Zehir, Jill P. Mesirov; Variant Review with the Integrative Genomics Viewer. Cancer Res 1 November 2017; 77 (21): e31–e34. https://doi.org/10.1158/0008-5472.CAN-17-0337
igv.js: an embeddable JavaScript implementation of the Integrative Genomics Viewer (IGV)
James T. Robinson, Helga Thorvaldsdóttir, Douglass Turner, Jill P. Mesirov
bioRxiv 2020.05.03.075499; doi: https://doi.org/10.1101/2020.05.03.075499

seeing some browser shots would help. In general, a few weird reads is no big deal. A lot of reads showing systematic issues might be more of a problem. It would also help to know what kind of reads. We trust Illumina more for base-level accuracy, and long reads more for contiguity. What species?

It's a human genome from the 1000 genomes project! And thank you, I just added a photo for clarification

thanks- can you share some of the track labels? I assume some of that stuff is gene tracks (e.g. green bars), but there are some coverage plots that don't show the read pileups, which makes me worry that I don't understand what's going on.

perhaps more to the point, there are 3 pertinent pieces of information you can get from this approach: (1) aligning reads against GRCh38, (2) reads against your assemblies, (3) your assemblies against GRCh38 (off the top of my head). I would suggest doing all 3 of these approaches, and comparing them. 1 and 3 can be visualized at the same time, which might be particularly helpful in finding e.g. regions of ambiguity due to repetitive regions. You appear to be doing 2, which might be a good way to identify regions to check by 1 and 3.

What do you mean by track labels? this site makes it seem like it's just the name of the file I'm using. In this case all my files are from de novo assemblies, I only got the FASTQ files off the internet

A filename might be helpful- basically I just don't know what I'm looking at. Even marking it up in MS paint or something would help.

Added! Thank you so much!


text:
Basic Docking Work Flow
I'm a computer scientist who is now getting into designing new molecules, with the goal that they dock well to a given target. While I'm used to coding up machine learning algorithms in python, I'm finding things a bit murky in terms of the how-to of docking. I lack pretty much any (medicinal) chemistry/biomedical engineering background and so understanding how to get things to work has been a real struggle.
In particular I've been looking at two githubs: this one using smina and this one using glide. Recently, my collaborators have steered me toward rDock/RxDoc.  Near as I can tell the basic ingredients for ANY docking software are:

a target file (rna or protein) in a .pdbdt,.pdb, or .mol2 file with
it docked to any ligand
a ligand file (mol2)
a target file with the ligand removed

But where I get lost is in any intuition on how to prepare these files, how to find them, and any intuition at all about how to make these work. I've seen suggestions of openbabel, etc, but there seem to be no tutorials online.
So, my question is: Suppose I have a molecule in rdkit and a protein or an RNA that I want to dock it to, accessed from the protein databank, say. What is the ``right'' way to do the following:

Download the receptor from pdb
Process the receptor to get file of the receptor with and without ligand
Process a molecule from rdkit so that the docking software can dock a (new) molecule to the receptor and give a docking score?

Ideally this should work for `any' receptor, but if the steps above need to be based on the specifics of the receptor, that would be helpful for me to know.
I'm shocked there isn't a beginners guide on this out there that's more easy to find, and I worry about going down this route without any intuition.

The basics of docking are explaining in several textbooks (eg. Essentials of computational chemistry). The science is rather important as there are lots of complications as protein are not solid blocks, water can form tight interactions (cf implicit water), Lipinski rules, the candidate compound needs to block the protein s activity (cf box size), Gibbs free energy (“score”) isn’t the sole value to sort on etc etc.

Regarding what to use, the Capri competition is a great way to see how programs compare. Glide scores well but it’s not open source. Autodock4 and Autodock Vina are very popular and highly documented —they are different in approach and the former is more adaptable and the latter is (misleadingly) simple. The papers describing these give a lot of info…

not sure I remembder it right, you put your target in a box , add the ligand somewhere in the box , minimize everything to find the most probable binding site ?, or you force the ligand in one particular site of the target ?


text:
Converting miRNA names to miRBase version IDs
I have a list of miRNAs IDs (2000-2500) that I want to find miRBase IDs for them. For example:
hsa-miR-106a ---> hsa-miR-106a-5p
hsa-miR-373* --> hsa-miR-373-5p
hsa-miR-33 --> hsa-miR-33a-5p
I used many R packages and websites such as (miRBaseConverter, miRNAmeConverter, MiEAA, miRSystem,miRandola) but non of them give me the comprehensive result.
e.g. For "hsa-mir-30c-1" miRNAmeConverter did not give miRBase ID but if you search it in miRBase it hsa two matures IDs: hsa-miR-30c-5p & hsa-miR-30c-1-3p.
miRBAseConverter is also like the mentioned one: It can not find the mature one for many miRNAs e.g. "hsa-mir-34" has to mature IDs: "hsa-miR-34a-5p" & "hsa-miR-34a-3p" but despite having them miRBaseConverter could not retrieve them.

Can anyone suggest any tools that might help to do this?


text:
Error in colSums(cts[[i]]) : 'x' must be an array of at least two dimensions
I'm getting an error when creating the counts from a single cell experiment.I'm identifying shared barcodes as well as lowUMI barcodes with the follow code:
shared = colnames(cts[[i]])[which(colnames(cts[[i]]) %in% lb[[libID]])]
lowUMI = colnames(cts[[i]])[colSums(cts[[i]])<1000]

but my code is erroring at
lowUMI = colnames(cts[[i]])[colSums(cts[[i]])<1000]

with this error
Error in colSums(cts[[i]]) : 'x' must be an array of at least two dimensions

My full code:
cts=list()
for (i in 1:length(mf$LibraryID)){
  # Read in all four input expression matrices
  message(i)
  message(mf$LibraryID[i])
  libID = mf$LibraryID[i]
  #lb[[libID]] = gsub("-1", "", lb[[libID]])
  lb[[libID]] = unique(lb[[libID]])
  ## by default we point to the 10X dir
  #data.dir = paste0(wd, mf$GSM[i], "/cellranger_out/",mf$LibraryID[i],"/filtered_feature_bc_matrix/")
  data.dir = paste0(wd2, mf$GSM[i], "/metadata/cellranger_ASI_run/",mf$LibraryID[i],"/filtered_feature_bc_matrix/")
      if (!dir.exists(data.dir)) data.dir = paste0(wd3, mf$GSM[i], "/", mf$LibraryID[i],"/outs/filtered_feature_bc_matrix/")
  if (!dir.exists(data.dir)) data.dir = paste0(wd2, mf$GSM[i], "/metadata/cellranger_ASI_run/",mf$LibraryID[i],"/outs/filtered_feature_bc_matrix/")
  if (dir.exists(data.dir)){
    message(data.dir)
    cts[[i]] <- Read10X(data.dir)
    #colnames(cts[[i]]) = gsub("-1", "", colnames(cts[[i]]))
    shared = colnames(cts[[i]])[which(colnames(cts[[i]]) %in% lb[[libID]])]
    lowUMI = colnames(cts[[i]])[colSums(cts[[i]])<1000]
    to_filter = intersect(shared, lowUMI)
    #to_filter = shared
    ok = colnames(cts[[i]])[!(colnames(cts[[i]]) %in% to_filter)]
    cts[[i]] = cts[[i]][, ok]
    colnames(cts[[i]])=paste0(mf$LibraryID[i], "_", colnames(cts[[i]]))
    if (i==1) {
      cts_all=cts[[1]]
    } else {
      cts_all = cbind(cts_all, cts[[i]])
    }
  }
}
dim(cts_all)

The error is clear - you are trying to use colSums on something with only one dimension, like a vector, whereas you need at least 2 dimensions, like a data.frame or matrix. But without any reproducible example, there is no chance that someone can fix the issue.

It does have 2 dimensions:
dim(cts_all)
106477  82029

Please can you provide a small sample of your data.

Might be wrong but it seems you are invoking colSums on the value of the colnames function which is a character vector so it has only one dimension. Are you trying to get the name of columns in cts where  colSums(cts[[i]])<1000 ?

@PaulEndymion Yes, I'm trying to retrieve all the barcodes that are less than 1000 and then filter out those barcodes.


text:
How do you replace multiple sequence names with new names using python?
I have into a dictionary some sequences. I would like to replace the sequences name for new names contained into a dataframe.
#my sequences into a dictionary
sequence = {"seq01": "agatcatggctaactatgcacgtgca",
           "seq02": "atatcatggctaactatgcacgtgca",
           "seq03": "tgatcatggctaactatgcacgtgca"}

#my data frame containing old names and new names
data = {"old_name": ["seq01", "seq02","seq03"], "new_name": ["sample1", "sample2", "sample3"]}
data = pd.DataFrame(data)
print(data)

I tried to replace the old names with a for loop, but it was imposible to do it.
Thank you

I am pretty sure there would be more elegant, in-place solutions but here is a quick solution to your problem. Should be OK if your sequence dictionary is not way too large.
seq_names = zip(data['old_name'], data['new_name']) # brings together old and new names in tuples
renamed_sequence = {}
for i in seq_names:
    renamed_sequence[i[1]] = sequence[i[0]]
renamed_sequence

EDIT: The section below is not quite accurate, please see the comments below.
And regarding your comment:

I tried to replace the old names with a for loop, but it was imposible
to do it.

You are not supposed to loop over and at the same time modify an element at the same time anyway (there might be side effects that you would not think of), hence a brand new dictionary in my code.

"You are not supposed to loop over and at the same time modify an element at the same time anyway": really? Why? Is that a Python thing?

@terdon, you are right, I should have been more explicit. For example changing the elements of a list while iterating over it would be fine, on the other hand, adding/removing elements while looping through a list would cause problems. In this particular case, the OP would be modifying the elements of a dictionary rather than changing its size, so an in-place solution using a for loop would be perfectly fine.

Ah, right. Things like entering an infinite loop with something like list.append(foo) running in a for i in list. Gotcha, thanks.

I would encourage you to improve your variable naming to something more descriptive: replace i with the tuple (old_name, new_name); names_tuple can become old2new or similar and new_dict could be called renamed. This yields the in my opinion much more readable assignment renamed[new_name] = sequence[old_name] which describes the code's intention. I think this is important because answers here serve as a reference for other coders and we should try to improve the overall quality of scientific software.


text:
how can I get expression of an inserted foreign genes?
hi we have transgenic mice with human gene inserted. if we profom rnaseq for the mice, how can we get the expression value of the human gene inserted to the mice? and use this expression to compare with other mice gene expression?
the mice ref genome doesn't have this piece of human gene sequence. anybody knows how?

You will need to create (and probably index) a "custom" genome/transcriptome, mouse + whatever human genes inserted, depending on the alignment method of choice. And how to do so would depend on the your tool of choice.


text:
How should I downsample and normalize R1s and R2s and incorporate this into Lexogen's QSPA tool
I am using Lexogen's Quantitative Sequencing Pool Analysis tool (here: https://github.com/Lexogen-Tools/quantseqpool_analysis) to analyze R1 and R2 files.
I have been able to successfully run this tool. However, I need to downsample and then normalize these results so I need to do this all again.
To downsample, I put the R1s (the ones generated after idemultiplexing) through the SeqTK_sample tool on UseGalaxy. Afterwards, I ran the QSPA tool with the new R1s and old R2s. I skipped the idemultiplexing part by # out the code. This is generating errors (screenshots are attached). I don't think these errors are from me # out the code but rather from the new R1s with the old R2s. I didn't receive any errors the first time I ran this tool.
I will be normalizing my data with DESeq2 using the summary_unique.tsv file.
If anyone has any tips/suggestions it would be much appreciated. I am sure there is a better way to do this. If you require anything from me just ask!
Cheers,

I don't know anything about this software in particular, but the traceback in the first screenshot makes it look like a compressed file is truncated.  QuantSeqPool calls cutadapt which in a roundabout way uses isal-l which uses isal.igzip which finally gives the "Compressed file ended" error.  I'd back up a step and look into what file is being given to cutadapt and why it's truncated.

Thanks for the advice, Jesse!
I compressed these files and loaded them onto the server. I then unzipped them using jar xf as conventional unzip methods weren't working. Do you think this could be an issue?

Sorry for the delay. I wonder if something went wrong during that process (truncated zip file?) and that might explain why "conventional" tools weren't working for you to unzip.  I'd suggest using a checksum program on both your local copy and the server copy (md5sum, sha256sum, whatever's available) and re-uploading any files that don't match.

I ended up figuring things out. I will post the solution for others to see.

The issue was that we downsampled with useGalaxy's SeqTK tool prior to running  UMI-tools and cutadapt (the 2nd and 3rd process, respectively, of Quantitative Sequencing Pool Analysis).
I tried downsampling with the R1 data, available after cutadapt (located in the "trimmed" folder), and everything ran fine. I believe the issues were caused by the barcode sequences being removed by the SeqTK tool prior to running UMI-tools.
You can (somewhat) customize the QSPA script by #'ing out parts of the script that you don't want to run. In this case, I #'d out Star Aligner, UMI-tools (deduplication), and Feature Counts so that the script would stop after cutadapt. After I downsampled with SeqTK, I ran the script using only the parts that I had previously #'d out.


text:
Help me in understanding the PDB file
Can anybody please explain the below line to me.
"We took the structure and coordinates of nogalamycin
from the X-ray structure determined in PDB code = 1D17"
What do I need to download from there....a pdb file or what?
Thank you

If you want to download the structure – yes, a PDB file or mmCIF file (different file formats with the same structure).

Thank you so much for your reply. But the pdb file is about the complex molecule(ligand-DNA). The above line says "We took the structure and coordinates of nogalamycin from the X-ray structure determined". The pdb id page shows ligand structure separately. .... download CCD file or instance file. Do you think I need to download only ligand file?

Ah, I overlooked that you are after the ligand only. Instance coordinates are taken from this structure (ligand extracted from the complex). Ideal coordinates are taken from the CCD (Chemical Component Dictionary) – the same ligand, but it's 3D structure doesn't need to be identical to the ligand in this complex. At least that's how I understand it. Another alternative is to download the PDB file with the whole structure and manually extract the ligand, but that's more cumbersome.

there is a https://models.rcsb.org/       ModelServer,
The ModelServer is a service for accessing subsets of macromolecular model data. I am sure there is a way to get just the ligand. Waiting for a Python client for the openAPI offered by it. somewhere it says Determines the output encoding (text based 'CIF' or binary 'BCIF'). Ligands can also be exported as 'SDF', 'MOL', or 'MOL2'

something like this https://github.com/williamgilpin/pypdb but for the modelserver

well openapigenerator docs states all operation on api with Return type Void so probably its not the right server ?

see https://bioinformatics.stackexchange.com/questions/8583/retrieve-id-ligand-from-pdb-file

Does this answer your question? Retrieve ID ligand from PDB file

Thank you everyone for your solution. I still did not figure out as I am totally new to this simulation field.

well wasnt able to figure out how to make openapigenerator work with rcsb pdb ModelServer https://models.rcsb.org/, I ended up with this code
copying the webserver manual input results for ligands:
import requests

def main():

  risposta = requests.get('https://models.rcsb.org/v1/1d17/ligand?encoding=cif&copy_all_categories=false&download=true&filename=out',

                          allow_redirects=True)

  open('out.cif', 'wb').write(risposta.content)

  print(risposta.status_code)

if __name__ == '__main__':

   main()

unfortunately the server gives you the wrong ligand (try it out)
it says :' /v1/{id}/ligand
Coordinates of the first group satisfying the given criteria.'
tried with:
risposta = requests.get('https://models.rcsb.org/v1/1d17/ligand?label_atom_id=NGM&encoding=cif&copy_all_categories=false&download=true&filename=out',

                       allow_redirects=True)

but no luck again, think problem is in downloaded file :
_model_server_params.name     atom_site
_model_server_params.value    '{"label_atom_id":"NGM"}'  ###### here

but not an expert of APis and openAPI at all
EDITED:
eventually I got it right using in the code above:
risposta = requests.get('https://models.rcsb.org/v1/1d17/ligand?label_comp_id=NGM&encoding=cif&copy_all_categories=false&download=true&filename=out',

                      allow_redirects=True)

so I get a out.cif file containing NGM
using another database and FTP :
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Nov 19 19:29:30 2021

@author: Pietro

https://bioinformatics.stackexchange.com/questions/18076/help-me-in-understanding-the-pdb-file#comment25289_18076

Ligand extract:

https://stackoverflow.com/questions/61390035/how-to-save-each-ligand-from-a-pdb-file-separately-with-bio-pdb

"""

from ftplib import FTP

import gzip

import os

from Bio.PDB import PDBParser, PDBIO, Select

def is_het(residue):
    res = residue.id[0]
    return res != " " and res != "W"

class ResidueSelect(Select):
    def __init__(self, chain, residue):
        self.chain = chain
        self.residue = residue

    def accept_chain(self, chain):
        return chain.id == self.chain.id

    def accept_residue(self, residue):
        """ Recognition of heteroatoms - Remove water molecules """
        return residue == self.residue and is_het(residue)

def extract_ligands(path):
    """ Extraction of the heteroatoms of .pdb files """

    for pfb_file in os.listdir(path):
        i = 1
        if pfb_file.endswith('.pdb') and not pfb_file.startswith("lig_"):
            pdb_code = pfb_file[:-4]
            pdb = PDBParser().get_structure(pdb_code, path + pfb_file)
            io = PDBIO()
            io.set_structure(pdb)
            for model in pdb:
                for chain in model:
                    for residue in chain:
                        if not is_het(residue):
                            continue
                        print(f"saving {chain} {residue}")
                        io.save(f"lig_{pdb_code}_{i}.pdb", ResidueSelect(chain, residue))
                        i += 1

def getFile(ftp, filename):
    try:
        # ftp.retrbinary("RETR", filename ,open(code+'.gz', 'wb').write)
        ftp.retrbinary("RETR " + filename ,open(code+'.gz', 'wb').write)
    except:
        print("Error somewhere !!!!")

code = '1D17'
filename = 'pdb'+code.lower()+'.ent.gz'                 ## pdb1D17.ent.gz

print(filename)

ftp = FTP('ftp.wwpdb.org')   # connect to host, default port

# ftp.login()               # user anonymous, passwd anonymous@
print(ftp.login())

ftp.cwd('/pub/pdb/data/structures/all/pdb/')

print(ftp.pwd()  )

getFile(ftp, filename)

zipped = gzip.open(code+'.gz', 'rb')

unzipped = zipped.read()

# print(unzipped)

saved = open(code+'.pdb', 'wb')

saved.write(unzipped)
saved.close()

path = './'

extract_ligands(path)

copying from : https://stackoverflow.com/questions/61390035/how-to-save-each-ligand-from-a-pdb-file-separately-with-bio-pdb
and using FTP from ftplib on http://ftp.wwpdb.org/pub/pdb/data/structures/all/pdb/   : PDB - FTP Archive over HTTP
NOGALAMYCIN is NGM third ligand downloaded (shuold be the only one according to https://www.rcsb.org/structure/1d17)
In any case I am failing at ubderstandung *,-)  how https://models.rcsb.org/
Parameters are related to https://pdb101.rcsb.org/learn/guide-to-understanding-pdb-data/beginner%E2%80%99s-guide-to-pdb-structures-and-the-pdbx-mmcif-format , the mmcif format
but as I told before I can't read inside openapi.json OpenAPI Specification


text:
Print first colum is the condition is true using awk
I have a binary file which looks like this,
Orthogroup  A   B   C   D   E   F
OG0000000   1   1   1   0   0   0
OG0000001   1   1   1   0   0   0
OG0000002   1   1   1   0   0   0
OG0000003   1   1   1   0   0   0
OG0000004   0   0   0   1   1   0
OG0000005   0   0   0   1   1   0
OG0000006   0   0   0   1   1   0
OG0000007   0   0   0   1   1   0
OG0000008   1   0   0   0   0   0
OG0000009   1   0   0   0   0   0
OG0000010   1   0   0   0   0   0

I want to use awk to print the first column if column D and E is 1 and rest is 0.
OG0000004
OG0000005
OG0000006
OG0000007

I tried to do this with individual files of each column  using the awk command
awk 'NR==FNR{arr[$0];next} $0 in arr' file1 file2

but doing this I had to generate multiple files and I thought using a binary file will be a more clean way.
Please Help,

It would be more helpful to paste the actual text instead of pictures of it.

I would do something like this:
$ cat test.tsv
groupA  1       1       1       0       0       0
groupB  1       1       0       0       0       0
groupC  1       0       1       0       0       0
groupD  1       1       1       0       0       0
groupE  0       0       0       1       1       1
groupF  0       0       0       1       1       0

$ awk '($5==1 && $6==1 && $2+$3+$4+$7==0){print $1}'
groupF

Thank you it worked

Straightforwardly
awk '$2==0 && $3==0 && $4==0 && $5==1 && $6==1 && $7==0 {print $1}' bin.file

A little spiffier
awk '$5 && $6 && !($2 || $3 || $4 || $7) {print $1}' bin.file

Thank you it worked


text:
Why does BWA MEM orientation contradict my library prep method
I have some RNA-seq data from a stranded paired end library prep, with dUTP and UDG preparation, so the orientation should be RF (confirmed with sequencing provider). I assembled the reads with Trinity, specifying RF for the strandedness, and I am mapping the reads to the transcripts with bwa-mem (I am aware bowtie2 or STAR might be better for this, but bwa is what I presently have installed, so just using it to check).
In the bwa output, I see lines # candidate unique pairs for (FF, FR, RF, RR): (0, 5372, 2, 0), which seems to contradict the expected orientation. Is there a reason that this is expected, or does it indicate a problem?

Not familiar with the protocol, but I guess RF means different things here. If I am right, with your RF protocol, the first read comes from the reverse strand of the transcript and the second read from the forward strand, but the pair orientation is still FR:
------------------->  transcript
  R2 ---->   <---- R1

The RF pair orientation should be
<----   ---->

You can confirm this by checking the mapping strand of read1.


text:
Remove gaps from alignment?
I have an MSA (protein sequence) and I have used various programmes (Clustal, Aliview, MEGA11 etc) to align. However, in all programmes I get many gaps which is not ideal as I am trying to construct a tree from this alignment.
My question is are there any tools to delete gaps in alignments?
Your help would be greatly appreciated.

You may also find useful the BMGE tool (stands for Block Mapping and Gathering with Entropy). BMGE selects regions in a multiple sequence alignment that are suitable for phylogenetic inference by computing their entropy-like score. This allows to distinguish among biologically expected and unexpected variability for each aligned character. BMGE also implements different types of trimming (for instance, deletes gaps) to minimize phylogeny reconstruction artifacts.
Available from here and also installable from Bioconda.

If you want an ungapped alignment, you could do an ungapped alignment (the older read aligners were ungapped, e.g. bowtie). However, I suggest that you don't want an ungapped alignment, you just want a stricter one, here's a few options...

Remove/trim poorly aligning sequences. You could try GUIDANCE, which has an explicit model or GBLOCKS for a more simple filtering criteria
Tune your alignment parameters for ClustalOmega/MAFFT etc. to make it stricter
Finally, consider that your sequences must be homologous otherwise you'll trim away sequences until you have a random alignment that appears good

This new tool, CIAlign, talks about some automated removal of gaps in multiple sequence alignments
https://twitter.com/thekatybrown1/status/1306215132670889984
https://github.com/KatyBrown/CIAlign

Gblocks is a tool for removing gaps/poorly conserved positions from alignments. Available here. Here is a webserver that runs it for you.
It's a bit old, possibly there are others at this point.


text:
How to extract part of protein complex separated with HEADER with two molecules from one PDB file
I have a PDB file which contain two molecules (receptor and ligand).
Each molecule will have its own header. All in ONE PDB file.
The header of receptor section looks like this (line 1-6 of the PDB file):
HEADER rec.pdb
REMARK original generated coordinate pdb file
ATOM      1  N   GLY A   1     -51.221 -13.970  37.091  1.00  0.00      RA0  N
ATOM      2  H   GLY A   1     -50.383 -13.584  37.482  1.00  0.00      RA0  H
ATOM      3  CA  GLY A   1     -50.902 -15.071  36.197  1.00  0.00      RA0  C
ATOM      4  C   GLY A   1     -49.525 -15.659  36.443  1.00  0.00      RA0  C

and this ligand section (line   11435 to 11440) of the PDB file
HEADER lig.000.00.pdb
ATOM      1  N   MET A   1      27.318 -26.957  12.663  1.00  0.00      LA0  N
ATOM      2  H   MET A   1      27.313 -27.570  11.870  1.00  0.00      LA0  H
ATOM      3  CA  MET A   1      28.374 -27.102  13.668  1.00  0.00      LA0  C
ATOM      4  CB  MET A   1      28.531 -28.564  14.090  1.00  0.00      LA0  C
ATOM      5  CG  MET A   1      27.224 -29.154  14.628  1.00  0.00      LA0  C

Note that the receptor and ligand section also contain the string RA0 and LA0
at the 11th column of PDB file.
What I want to do is to rename the chain of receptor as chain A and ligand as chain B.
In order to do that I intended to extract both part into two different objects first,
then rename the chain and finally put them together again.
I tried this code with Bio3D.
But doesn't work:
library(bio3d)
pdb_infile <- "myfile.pdb"

pdb <- read.pdb(pdb_infile)
receptor_segment.sele <- atom.select(pdb, segid = "RA0",  verbose = TRUE)
receptor_pdb <- trim.pdb(pdb, receptor_segment.sele)
ligand_segment.sele <- atom.select(pdb, segid = "LA0",  verbose = TRUE)
ligand_pdb <- trim.pdb(pdb, ligand_segment.sele) # showed no entry

What's the way to do it?
I'm open to solution with R or Python.

Since this is not a correct PDB file, it can be interpreted differently by different PDB reading libraries. The most reliable way to change it would be to read it as plain text and change the chain ID to B in the lines that have segment ID LA0.
Note: Strings RA0 and LA0 start at column 73. Columns in the PDB format correspond to characters. This is not a whitespace-delimited file.
Chain ID is in column 22.
Segment ID is in columns 73-76.
As a side note, segment ID was only briefly part of the official PDB spec (around 1996), nevertheless it's supported by many programs.
Of course, if you have only one file that you want to modify, the fastest way is to do it manually in a text editor.


text:
Generating intercalation site in DNA for custom sequence
Thank you for your help.
Can anybody please tell me how I can generate an intercalation site in DNA at the base pairs I am interested on. I appreciate if you can please guide me to any tutorials or videos which shows how to do that. Actually, I want to intercalate at particular site and I did Autodock to but it shows that the ligand is not at all intercalated.

As mentioned in a previous Q intercalation requires pi stacking and interactions of quadrupole moments. The Lenard-Jones potential is a poor proxy for this and you may need to use quantum mechanics —have you looked into Gaussian etc.? (I may be wrong though if there are lots of polar interactions)

Link to previous question which this may be a duplicate of: building complex drug-dna for AMBER software


text:
HG19 Position meaning
I am looking at genomic data (HG19). I have many SNPs, their chromosomes and positions.
I want to look at certain SNP (suppose its chromosome is 1 and position 77,226,919), and extract all SNPs in distance 1mb from that SNP.
Are the HG19 position in terms of base-pairs? Meaning that I need to take all SNPs in position $77,226,919 \pm100,000$? Or is it something else?

A megabase is 1_000_000; not 100_000.

Yes, they are base pair positions

HG19 refers to the build used as the reference genome. To extend on what @Emily_Ensembl said, base pair positions may change between different builds.
Please look at the following link for a more detailed explanation: https://gatk.broadinstitute.org/hc/en-us/articles/360035890951-Human-genome-reference-builds-GRCh38-or-hg38-b37-hg19


text:
What does 'Human Alternative sequence Gene' mean in Ensembl?
I have downloaded a dataset containing RNA-Seq RPKM (reads per kilobase per million) values for 52376 genes (in humans). To find RPKM values for the gene KIR2DL2 in the dataset, I have to search for it using its Ensembl ID.
After searching for this gene using the Ensembl browser, I obtain many results.

For all the results in the first page, it says Human Alternative sequence Gene after the gene name. When I have searched for other genes using their Ensembl ID I get only one result:

What does Human Alternative sequence Gene mean in Ensembl? I am not sure how to tell which is the correct Ensembl ID for the KIR2DL2 gene.
Thanks.

These are genes on the patches and haplotypes. Patches are repairs in the genome, where they have discovered the sequence or assembly is incorrect, so there is an alternative sequence that is created on top of the reference assembly (the reason they don't change the reference is that then there's a new reference and this causes a massive headache for everyone). Haplotypes are alternative versions of genomic loci that exist in some individuals, such as the MHC locus, which are different and common enough to require a proper assembly change, more than just genetic variants.
The genes that are found on the alt sequences are listed as Human Alternative sequence Gene. They are both real genes and correct for the gene name.
If the gene is on a patch (the location will include the word PATCH in its name) you should probably use the patch gene because the primary assembly gene may be dodgy.
It gets more complicated if you're working with haplotypes, as all the versions are correct. Most people take the easy option of using just the primary assembly, but if you're variant calling, it may be worth trying to determine which haplotype the individual has and calling against that. This is particularly true for highly variable regions like MHC, where there are genes only found on the alternative haplotypes and not found on the primary assembly.
There's more information about alternative sequences from Ensembl here. And here from the GRC.


text:
replacing periods with hyphen
I am a software engineer and not a bioinformatician. I am looking at some code where they are looking at aligned proteins that are in the A2M alignment format (https://compbio.soe.ucsc.edu/a2m-desc.html).
According to the page, the
Uppercase characters and "-" represent alignment columns,
and there must be exactly the same number of alignment columns in each sequence

Now, the code I am looking at is doing something like:
sequence = sequence.replace(".","-")

So replacing the dots with hyphen and I am not sure what the purpose of that might be.
For anyone interested, the code is here (See line 182): https://github.com/debbiemarkslab/DeepSequence/blob/master/DeepSequence/helper.py

In A2M format, upper case letters represent matches, lower case letters represent inserts, dashes represent deletions, and dots (or spaces outside the identifier lines) represent gaps aligned to inserts.
So, both "-" and "." are essentially gaps, but assumed to have different origin. This information is supplementary, and most MSA programs ignore dots and spaces or treat them incorrectly.
The code you're looking at, in lines 187–197, removes sequences that contain characters not from self.alphabet (defined in lines 74–86) by checking if letter not in alphabet_set and letter != "-" (line 192). Therefore, just for the sake of parsimony, it looks natural to transform "." into "-" beforehand.

Thanks for the explanation!


text:
Trinity Multi-species compatibility
I am using Trinity from last couple-up months or so. Initially, using same species in different condition was worked fine for me. Now i am trying to do a DE analysis using multiple species. Actually, i have 3 fish performs similar behaviour  and another 3 performs similar strategy. Here i considered 3 similar species as replicate for one type (and totally i have two type) . Though these are 6 different species but belong to same family.
I have already made a combine assembly according to trinity protocol (--samples_file). After assembly i have checked the separate alignment rate for each species and it was above 85%. And lastly i checked heat-map and it gave me a result that somehow match with my expectation.
But i am a bit confused about the possible generation of "assembly chimera". In most of the cases i found people do DE analysis with same species. Also didn't find any reference article that did multiple species DE analysis using Trinity pipeline. So it will be a great help for me if i get some suggestion that either i should go with existing result or i need to find another way.

Given that you're using multiple species / samples, the results from this are going to be quite weird.
For doing such a comparison, I think the best approach is to do a combined assembly (as it looks like you've done) to generate a single chimeric transcriptome (which will not be a good representation of any of the component species), then do a differential expression analysis using that combined assembly as the reference.
Ideally you would have multiple replicates for each different species, and then add the species and conditions into the statistical model, so that any species-specific variation can be identified and/or filtered out.
Six separate assemblies that are later merged will be less useful for a differential expression, because it will bias towards differences in the assemblies. This is a similar rationale to that used for generating UMAPs from multiple single-cell datasets where a combined UMAP is generated for the full dataset with identical parameters, so that sample-specific differences in the UMAP parameters are less of an issue.
I think it will be okay to go ahead with the combined assembly, bearing in mind that the lack of replicates for the differential expression study is going to make it more difficult to identify biological significance. Using split assemblies will not make that better; if anything, it'll make it worse.

Thank you very much.


text:
At what point is a gene different between species that we call it a different gene?
I am new to genetics and I know humans share many genes with mice for instance but that there are slight differences in conserved nucleotide sequences.  Is there a community consensus around at what point a gene in one species is a different gene?  How are regulatory elements like enhancers factored into this definition?

From an evolutionary angle, you are describing an ortholog: "Orthologs are genes in different species that evolved from a common ancestral gene by speciation, and, in general, orthologs retain the same function during the course of evolution." quoting here. A lot of effort has been made to construct ortholog databases (e.g. see quest for orthologs). They are distinguished from paralogs (approximating duplicated genes) and lateral gene transfer events, there is however, no magic cutoffs for defining orthology.
Your 2nd point about regulatory elements, well they could fit in the same framework, any sequence can be clustered. What is a gene anyway? Is a noncoding RNA a gene? How do we go from a collection of transcripts to this unit we term gene? There is some ambiguity...


text:
Visualizing top expressed genes
I have a spreadsheet of CPM values
I want to visualize the top 20 genes expressed across my samples
This is how my data look like
> head(b[1:5,1:5])
                 30566-001 - CPM 30566-002 - CPM 30566-003 - CPM 30566-004 - CPM 30566-005 - CPM
hsa-miR-548ap-3p               0               0               0               0               0
hsa-miR-548t-3p                0               0               0               0               0
hsa-miR-548e-3p                0               0               0               0               0
hsa-miR-548az-3p               0               0               0               0               0
hsa-miR-548f-3p                0               0               0               0               0
> dim(b)
[1] 2632   42
>

How I can do that?

Just make a heatmap?

Then from a data frame with 2632 genes, how I select a list of genes for heatmap? Because a heatmap with 2632 genes is not clear while my main purpose is to shows a list of highly expressed genes across my samples

select the N genes with the highest average expression across all samples?

Gene selection is entirely on you and on the question you ask. Most expressed genes would be selected by rowMeans or rowMedians or even rowSums but this is not very informative as digital expression levels (=counts) is influenced by many technical factors and there is (to my knowledge) no direct link to biological relevance.
Alternatively, you can select by genes that vary a lot between experimental groups. Naively that could be the rowwise variance (rowVars) or standard deviation of the logCPMs (use logscale, otherwise variance is a function of the expression level) or the output of a differential analysis. Ideally plug the raw counts into DESeq2 or similar frameworks. If you only have CPMs then you probably want to feed the logCPMs into limma-trend pipeline, see for a discussion on ways to get going with normalized counts for differential expression in this Bioconductor thread.
Visualization is then again on you. Either make a heatmap (see for example ComplexHeatmap package) or aggregate the genes into a PCA (see for example PCAtools package at Bioconductor), all mentioned packages contain extensive example code and background information which is fully sufficient to run such an analysis.


text:
How to split genotypes into two groups from a SNP txt file contains -1, 0 and 1 in R?
I have a SNP file(it is 96 obs. of 1178 variables with -1, 1 and 0 and in some places NA) and I need to calculate Jaccard genetic distance from them in R. Is there any way to convert these 3 values to 0 and 1? I tried to change from -1 to 2 before for a different distance, but the Jaccard distance is binary and you cannot assign 2. Is there a way to re-code Jaccard to take into account -1 and give a good distance matrix? (Jaccard must remain binary and contain only 0 and 1). I will need it later to divide the subpopulation into 2 groups.
Here is an example input data file:
https://pbgworks.org/sites/pbgworks.org/files/snp.txt

THis should be pretty simple to achieve, but could you clarify exactly your rules - i.e. which number gets converted to which?


text:
Where can I get info about A. thaliana proteins' locations?
I have a list of A. thaliana proteins and I'd like to know their location in the cell. I've found SUBA but they seem not to share their dataset and only allow queries through their website. It's inconvenient for me, because I need to locate 1.6k of proteins and process the results further. Is there a dataset or an API I could use instead?

When cross-posting on multiple sites, it is considered polite to include links to the other posts as a courtesy to members (e. g. https://biology.stackexchange.com/questions/105381/where-can-i-get-info-about-a-thaliana-proteins-locations)

Thanks for the advice, I'll do it next time.

Given the fact that it is a simple query, in my opinion, the best approach would be an advanced search on AmigoGO. There is no need to bother with complex methods like programmatic access APIs, as you basically have only two criteria to meet: organism and cellular localization.
After you specify those, you can download the file using the appropriate button, or mess a bit with the columns that you are interested in.
Once you have the file, you can then manipulate it further to your specific need with any tool of your preference, for example, to retain only proteins of your interest.
Hope this helps.


text:
Publicly available genotype / phenotype dataset?
I'm looking for publicly available genotype data (eg. 1000 Genomes or HapMap) that also has associated phenotypes (any traits).  I'm wanting to use it as a test run for a GWAS pipeline. Does anyone know if this exists anywhere?

See openSNP here:
https://opensnp.org/
It is a web site for people to upload their own genetic datasets, providing freely-accessible public data with annotated phenotypes.

openSNP allows customers of direct-to-customer genetic tests to publish their test results, find others with similar genetic variations, learn more about their results by getting the latest primary literature on their variations, and help scientists find new associations.

The large collaborative research genetic data projects typically have very strict access constraints due to the realisation that genetic data is personally-identifiable data, and the expectation of anonymity. That means there's a low likelihood that a large research project will make their phenotype + genotype data publicly available.


text:
Changing the order of colors in pheatmap
I plotted a heatmap
I want cancer to be red and healthy to be blue
   pheatmap(heat1, annotation_col=df,color=colorRampPalette(c("navy", "white", "red"))(50))

You are seeing cancer is in blue in annotation, how I can change that

It is possible with annotation_colors argument :
First create a list with your conditions and then add the argument :
annot_colors=list(HealthStatus=c(Cancer="#F0978D",Healthy="#63D7DE"))

pheatmap(heat1,
         annotation_col=df,
         color=colorRampPalette(c("navy", "white", "red"))(50),
         annotation_colors=annot_colors)


text:
How to identify in a RNA seq data in which sample is in which cell using R using qusage
I am using RNA seq data and have been using QuSage in R in hopes to try to identify different samples and figure out which cell it belongs to.
I am trying to visualize and utilize the results but I am having trouble doing that. Does anyone have an idea?

QuSAGE
Quantitative Set Analysis for Gene Expression method
https://bioconductor.org/packages/release/bioc/html/qusage.html

Please be more specific. What have you tried, and where are you running into problems?


text:
Updating the GFF3 + Fasta to GeneBank code
I'm trying to convert gff3 and fasta into a gbk file for usage in Mauve. I've found a solution but the code is outdated:
"""Convert a GFF and associated FASTA file into GenBank format.

Usage:
gff_to_genbank.py <GFF annotation file> <FASTA sequence file>
"""
import sys
import os

from Bio import SeqIO
from Bio.Alphabet import generic_dna
from BCBio import GFF

def main(gff_file, fasta_file):
    out_file = "%s.gb" % os.path.splitext(gff_file)[0]
    fasta_input = SeqIO.to_dict(SeqIO.parse(fasta_file, "fasta", generic_dna))
    gff_iter = GFF.parse(gff_file, fasta_input)
    SeqIO.write(gff_iter, out_file, "genbank")

if __name__ == "__main__":
    main(*sys.argv[1:])

But the Bio.Alphabet module is removed now. I tried to read the documentation on these changes and how to update the code (https://biopython.org/wiki/Alphabet) but I'm still confused about how to update the code. If I remove the Bio.Alphabet and the generic_dna, the code will return ValueError
ValueError: missing molecule_type in annotations

I know that I have to type the molecule_type = "DNA" but I'm confused as to where should I put it because the code directly inputs the fasta_file into the SeqIO.parse() function.
Anyone knows how to update code related to SeqIO.parse() with the removal of Bio.Alphabet? Or anyone knows how to convert GFF3 + Fasta into gbk files?
Note: I can't use the seqret because my files is just too big for it, same problem as this question.
Thanks in advance!

You can set the molecule type on the existing annotations attribute of the record objects just before they get to SeqIO.write.
There might be an even leaner way to do it, but if open the output file (so that you can do multiple writes to one file handle) and then iterate over each record from the GFF iterator, you can do it like this:
with open(out_file, "wt") as f_out:
    for record in GFF.parse(gff_file, fasta_input):
        record.annotations["molecule_type"] = "DNA"
        SeqIO.write(record, f_out, "genbank")

(The documentation is vague on that point but if you check you'll see that record.annotations is a dictionary that you can modify to set the molecule type.)

not sure it's leaner but something like SeqIO.write((SeqRecord(r,annotations={"molecule_type": "DNA"}) for r in GFF.parse(gff_file, fasta_input)), out_file, "genbank")

Thanks Jesse, it works like charm!


text:
bowtie No alignments
This question was also asked on Biostars
I am working on miRNA alignment using bowtie. But it seems to fail to map. Here are the commands:
bowtie-build -f miRNA.fa miRNA
bowtie --norc -k 1 -v 0 /public/home/wschen_gibh/cpa-seq/ref/miRNA/miRNA -1 SRR16643593_1.NEW.trim.fq -2 SRR16643593_2.NEW.trim.fq -S test.sam

and here is the result:

and the output sam file:

Ah, this is "CAP-seq", so there are presumably UMIs in those sequences that you need to deal with first

it would help to see the reference file. Can you confirm that you actually expect anything to align to that reference?

You should be using Bowtie 2 instead of bowtie:

Bowtie 1 was released in 2009 and was geared toward aligning the relatively short sequencing reads (up to 25-50 nucleotides) prevalent at the time. Since then, technology has improved both sequencing throughput (more nucleotides produced per sequencer per day) and read length (more nucleotides per read).

http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml#how-is-bowtie-2-different-from-bowtie-1
In particular, Bowtie 2 supports local alignment, which is essential when mapping microRNA. All microRNA sequences are short, which means that there will be adapter read-through in the sequenced reads. This adapter sequence will not map to the genome, and that will cause global aligners (like Bowtie) to fail.
An alternative approach would be to trim off palindromic adapter sequences prior to aligning (which it is possible you've tried to do given the file name). Because of adapter read through, standard adapter trimming will not be sufficient. What needs to be trimmed off is the reverse-complement of the adapter sequences. I have no experience using fastp, so don't know if it is capable of doing that.
However, even after trimming there will still be some additional short sequences left on some reads, so a local alignment would still be a good idea.
Doing a local alignment with Bowtie 2 is much easier and faster than trying to troubleshoot adapter read-through issues. The command line you would run is similar, but not identical:
bowtie2-build -f miRNA.fa miRNA
bowtie2 --local -x miRNA -1 SRR16643593_1.NEW.trim.fq -2 SRR16643593_2.NEW.trim.fq | samtools sort > test.bam


text:
Phasing genotyped data using shapeit: ERROR: Wrong ordering in physical positions
I am to phase my genotyped data with Shapeit (from ped or bed files) and I keep on getting this error message:

ERROR: Wrong ordering in physical positions curr_pos=13915
prev=157299737

Any idea on how to fix this?

From the website:
"SHAPEIT is a fast and accurate method for estimation of haplotypes (aka phasing) from genotype or sequencing data. "
https://mathgen.stats.ox.ac.uk/genetics_software/shapeit/shapeit.html

It appears that your SNPs aren't ordered correctly - can you post some of those lines from the .bim file? Also, it seems like you are using shapeit2. Unless you have a good reason, you should be using the most up to date version, shapeit4.


text:
Definition of Copy Number Abundance
Reading from this paper:

The remaining 9 new copy-number alterations showed altered coverage patterns on whole-genome sequencing but could not be confirmed by orthogonal methods because of their small size, low abundance, or both (Fig. 2C, Fig. S2C, and Table S5)

In the supplementary file, CNA abundance as percentage is defined as:
$$Abundance = \frac{2^{L2R} - 1}{\frac{CN}{2} - 1}*100$$
where L2R is the log$_{2}$ normalized coverage ratio vs. a panel of normals and CN is the estimated copy number for the event.
I haven't found a better definition/explanation behind "CNA abundance", is it defined elsewhere, or is it a common term in CNA analysis?

My guess by analyzing the definition is that CNA abundance is a proxy for the percentage of (tumor) cells carrying the alteration.


text:
Polishing PacBio or ONT with Illumina
Which tool would be good to polish PacBio or ONT with Illumina?
Thank you in advance

Do you mean to correct PacBio/ONT reads or polish their assembly?

Metagenomic or single genome?

PacBio and ONT are very different technologies; it would be helpful to be more specific. Regardless, I can only give my understanding as it relates to ONT reads....
TLDR: use medaka. Don't use Illumina reads.
Nanopore reads called using the most recent super-accuracy basecallers have a very high consensus accuracy (q40-q50) for SNPs and non-homopolymeric INDELs after nanopore-only correction. ONT provide their own tool, medaka, which uses assembly-mapped nanopore reads to correct the assembly. Medaka uses the base calling model from Guppy to predict the most likely sequence, given the input reads and the assembly.
If you're using Illumina reads to polish a nanopore assembly, it should only be used to correct INDEL errors (and ideally, only homopolymer errors). The two most commonly-used tools that can correct nanopore reads using Illumina reads are Racon and Pilon. Based on the recency of the commit histories, I'd guess that Racon would be the most likely to work well. It looks like Racon can do both assembly correction and read correction.

Why don’t you use a hybrid assembly which pretty much does it by default? Tools like OPERA-MS assemble the short reads with MEGAHIT and then uses nanopore to build bigger contgs. Otherwise, I recommend pilon which probably the most popular tool.

For PacBio HiFi reads (or other techs of comparable quality), it is not recommended to polish with Illumina at all. You run the risk of collapsing true heterozygosity/variation in your underlying genome and yielding a haploid assembly that does not represent either haplotype.
If that underlying variation is an issue because you need a haploid assembly, you are free to use a deduplication technique such as purge_haplotigs or similar (reviewed here) to obtain a primary set of contigs/scaffolds. Assigning p and q haplotypes is something that you can do with a tool like falcon-phase (if you have Hi-C) or a similar technique (I haven't been following the literature for a couple years, possibly there's something a lot better now!).
For a recent review of recommended strategies you can read the T2T consortium benchmarking paper on the topic.