text:
curl: (22) The requested URL returned error: 400 Bad Request when using entrez direct
I've just dowloaded entrez direct on my Mac and I was trying to perform some simple queries like:
esearch -db pubmed -query "tn3 transposition immunity" | efetch -format medline

but I keep getting this error
curl: (22) The requested URL returned error: 400 Bad Request
 ERROR:  curl command failed ( Mar 19 Ott 2021 13:02:26 CEST ) with: 22
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi -d db=pubmed&id=Unable%2Cto%2Clocate%2Cxtract%2Cexecutable.%2CPlease%2Cexecute%2Cthe%2Cfollowing%2Cnquire%2Cdwn%2Cftp.ncbi.nlm.nih.gov%2Centrez%2Centrezdirect%2Cxtract.Darwin.gz%2Cgunzip%2Cf%2Cxtract.Darwin.gz%2Cchmod%2Cx%2Cxtract.Darwin&rettype=medline&retmode=text&tool=edirect&edirect=16.1&edirect_os=Darwin&email=mariocatalano%40MBP-di-Mario.lan"

I'm using bash scripts on my Mac's terminal
How should I fix this?

update your curl or esearch version, use the miniconda/brew versions. If you already have conda, you can create an environment and install it on it.

Looks like you haven't installed it correctly. The error is saying it cannot locate/extract the executable.

@zorbax -- it is recommended to use the official version of Edirect instead of the  ones from miniconda, brew, etc as they often tend to be quite behind the latest, official version.

It's not a curl problem.  (The 22 error is just passing along the server-side HTTP 400 error, which is because it dumped an error message into a URL.)  I think terdon is right and there's something wrong with your edirect installation.  You have esearch and efetch on your PATH but apparently not xtract and/or xtract.Darwin.

More detail:
Look toward the end of the error message, at the text being passed to curl's -d argument:
db=pubmed&id=Unable%2Cto%2Clocate%2Cxtract%2Cexecutable.%2CPlease%2Cexecute%2Cthe%2Cfollowing%2Cnquire%2Cdwn%2Cftp.ncbi.nlm.nih.gov%2Centrez%2Centrezdirect%2Cxtract.Darwin.gz%2Cgunzip%2Cf%2Cxtract.Darwin.gz%2Cchmod%2Cx%2Cxtract.Darwin&rettype=medline&retmode=text&tool=edirect&edirect=16.1&edirect_os=Darwin&email=mariocatalano%40MBP-di-Mario.lan"

It's putting this in the id= part:
Unable%2Cto%2Clocate%2Cxtract%2Cexecutable.%2CPlease%2Cexecute%2Cthe%2Cfollowing%2Cnquire%2Cdwn%2Cftp.ncbi.nlm.nih.gov%2Centrez%2Centrezdirect%2Cxtract.Darwin.gz%2Cgunzip%2Cf%2Cxtract.Darwin.gz%2Cchmod%2Cx%2Cxtract.Darwin

which is maybe something like (trying to guess which characters got turned into %2C):

Unable to locate xtract executable. Please execute the following: nquire -dwn ftp.ncbi.nlm.nih.gov entrez/entrezdirect xtract.Darwin.gz; gunzip -f xtract.Darwin.gz; chmod +x xtract.Darwin

But you should already have xtract and xtract.Darwin from the same edirect install that gave you esearch and efetch.
For what it's worth the official Entrez Direct install on Mac OS looks a little goofy.  When I try it myself it ends with:
Entrez Direct has been successfully downloaded and installed.

-e In order to complete the configuration process, please execute the following:

  echo "export PATH=\\\${PATH}:" >> $HOME/.bashrc
  echo "export PATH=\\\${PATH}:" >> $HOME/.zshrc

or manually edit the PATH variable assignments in your .bashrc and .zshrc files.

Would you like to do that automatically now? [y/N]

Those commands don't look right (what's getting added to the path? an empty string?) and if I say y, I get this at the end of .bashrc and .zshrc which definitely isn't right:
export PATH=${PATH}:

It seems like an empty string at the end of PATH behaves like . (current directory), so I can go into ~/edirect and run your command and it does actually work.  So this isn't exactly your problem... but still, that seems wrong.

Note that macOS uses login shells by default for historical reasons, which means that if the OP is using bash (which is no longer the default on macOS), then ~/.bashrc is irrelevant and they should use ~/.profile instead. Either way though, that command is quite wrong, as you point out.

This was an issue with the Edirect installation script. The install-edirect.sh script has been updated to fix this issue. Please follow the instructions here to install.

Thanks!  The install still seems to misbehave on Mac OS, though.  When I say ymy ~/.bashrc is now untouched and my ~/.zshrc contains literallyexport PATH=${PATH}:/Users/jesse echo export PATH=${PATH}:/Users/jesse

@Jesse thank you for trying it out. This is yet another issue that needed a quick fix. The updated script is already on FTP so if you try again it should work. Sorry for the inconvenience.

No trouble, thanks for working on this.  Still behaving strangely for me but when I look at the script I think I see one reason.  It references DIR when setting up PATH, but that's the original working directory / the directory containing the install script, not the destination directory, right?  Should that instead be a $PWD or something once inside the ~/edirect location, maybe? We might be drifting from OP's question but I'm happy to help debug this via email or a github gist or something.

@Jesse, turns out this is yet another issue that had to be fixed. The new installation script was derived from the old setup.sh file, which was run within the edirect folder. In that context the DIR variable was correct. The new code now sets DIR after calling cd edirect, so its value is now also correct. Sorry for the inconvenience, and thanks again for your reports.

Yep, works as expected now on Mac OS.  Thanks!


text:
How can I calculate the distance w.r.t. a beta sheet?
Suppose, I have a protein (say, 16PK).
and, I want to calculate

the distance between one of its alpha-carbons and the nearest beta-sheet.
the distance between two strands in a beta-sheet (i, j).

How should I proceed to achieve this?

The documentation for the PDB module in biopython, which you have tagged, is very clear about how distances are measured between two given Cα atoms.
To find out the secondary structure and interactions between strands of a sheet, you may want to look at the header of the PDB. This is the stuff before the ATOM/HETATM entries. SHEET entries are what you are after. Please read the PDB format specifications for secondary structure. Biopython PDB does not parse headers well, but reading the file and finding the SHEET entry and positionally splicing out the strings of interest from the entry is beyond straightforward.
Example take from the above link:
         1         2         3         4         5         6         7         8
12345678901234567890123456789012345678901234567890123456789012345678901234567890
SHEET    1   A 5 THR A 107  ARG A 110  0
SHEET    2   A 5 ILE A  96  THR A  99 -1  N  LYS A  98   O  THR A 107
SHEET    3   A 5 ARG A  87  SER A  91 -1  N  LEU A  89   O  TYR A  97
SHEET    4   A 5 TRP A  71  ASP A  75 -1  N  ALA A  74   O  ILE A  88
SHEET    5   A 5 GLY A  52  PHE A  56 -1  N  PHE A  56   O  TRP A  71
SHEET    1   B 5 THR B 107  ARG B 110  0
SHEET    2   B 5 ILE B  96  THR B  99 -1  N  LYS B  98   O  THR B 107
SHEET    3   B 5 ARG B  87  SER B  91 -1  N  LEU B  89   O  TYR B  97
SHEET    4   B 5 TRP B  71  ASP B  75 -1  N  ALA B  74   O  ILE B  88
SHEET    5   B 5 GLY B  52  ILE B  55 -1  N  ASP B  54   O  GLU B  73

what does it mean by i and j in case of beta-sheets?

It likely depends on the context of the paper or textbook... But in the PDB format snippet above sheet A is composed of 5 strands, so the mean distance between the closest Cα atoms in strand i and those in j varies depending on their adjacency...

source: Stack Exchange

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

sort: votes


text:
When using cutadapt, can I specify an R2 adapter as optional when I have a required R1 adapter?
I'm using cutadapt 3.5 to trim adapters and perform some filtering on paired-end data.
Both R1 and R2 sequences have 3' adapters that might be found depending on the sequence length, but R1 also has an anchored 5' adapter that I want to require for the pair to be kept in the output.  For example, if that required R1 5' adapter is AAAAAAAA, the R1 3' adapter is TTTTTTTT, and the R2 3' adapter is GGGGGGGG, and I have this input:
R1:
@read1 should be kept
AAAAAAAACGTCCTGGGATTTGTAATAATATTTTAGTTCTGAGCGACAAGTAAGGGATAATTTTTTTT
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@read2 should be filtered
TACCAACTACATTTAGCTTCAGGCTAGTGATGCCCGCCGTCGGCACACTGGACACATGGCTTTTTTTT
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@read3 should be kept
AAAAAAAACAATTTTTACTCTAGAAATGTCTGTGCTCATCACCCCGACACCGAATAGCTATGACTCAC
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

and R2:
@read1
AGGCTCACCGCCACTGTTGTACCTTCCTATCGCCACTCTAAGATCTATGTAACCTCTCCCGGGGGGGG
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@read2
TGAATAATACTGCGGTTTAATCGATATATTCGGGATTATGCAAGACACCCTACGTACTTAGGGGGGGG
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@read3
TACTGGAGGGTACAAGCCCGGACTATCCACGGCGTCAGCTGGCTTAGCATTGAAGGTCGGCGGTGTGT
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

I'm aiming to get this R1 out:
@read1 should be kept
CGTCCTGGGATTTGTAATAATATTTTAGTTCTGAGCGACAAGTAAGGGATAA
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@read3 should be kept
CAATTTTTACTCTAGAAATGTCTGTGCTCATCACCCCGACACCGAATAGCTATGACTCAC
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

And this R2 out:
@read1
AGGCTCACCGCCACTGTTGTACCTTCCTATCGCCACTCTAAGATCTATGTAACCTCTCCC
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@read3
TACTGGAGGGTACAAGCCCGGACTATCCACGGCGTCAGCTGGCTTAGCATTGAAGGTCGGCGGTGTGT
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

This approach doesn't work as it discards read3:
cutadapt -a '^AAAAAAAA...TTTTTTTT' -A 'GGGGGGGG' example.R{1,2}.fastq --discard-untrimmed -o out.R1.fastq -p out.R2.fastq

I expect that's because the required/optional behavior only applies to the linked adapter for R1, and R2 is considered untrimmed for read3, which then filters the whole pair because of the default --pair-filter=any behavior (if either read is untrimmed, the pair is discarded).  I found I can get the filtering I want with --pair-filter first, but then the other filtering criteria like quality and length will only apply to R1, too.
My question is, can I get the trimmed/untrimmed behavior I'm looking for (read1 and read3 trimmed, read2 untrimmed) while also applying cutadapt's filtering to both reads in each pair? In other words, can I specify to cutadapt that R1's anchored adapter is required and the other two optional, without resorting to changing the pair filter setting?

This is not possible in a single cutadapt call, but can be done fairly seamlessly with a pipe between two calls and with cutadapt's interleaving feature.  (Thanks to the author for help with this.)

Giving the R2 3' adapter to one cutadapt command and interleaving the read pairs to piped output, a second cutadapt command can take that interleaved output and use both of R1's adapters as a linked adapter:
cutadapt -A GGGGGGGG --interleaved example.R{1,2}.fastq | \
cutadapt -a "^AAAAAAAA...TTTTTTTT" --interleaved --discard-untrimmed -o out.R1.fastq -p out.R2.fastq -

Now I can add other filters to the second call.
An additional read pair that should be filtered:
R1:
@read4 would be kept but R2 is too short
AAAAAAAACAATTTTTACTCTAGAAATGTCTGTGCTCATCACCCCGACACCGAATAGCTATGACTCAC
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

R2:
@read4
CTACTAGACATATACCCTGGGGGGGGTTCCCCGTGAAGACGGTTACGATGTGGACAAGCACTCCAGCG
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

I can now give a minimum length that applies to both in each pair:
cutadapt -A GGGGGGGG --interleaved example.R{1,2}.fastq | \
cutadapt -a "^AAAAAAAA...TTTTTTTT" --interleaved --discard-untrimmed --minimum-length 30 -o out.R1.fastq -p out.R2.fastq -

With that I get read1 and read3 in the output, as intended.

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


text:
How to Find List of Genomic Coordinates from List of Genes using UCSC Website?
Suppose I have a list of genes and I want to find their genomic coordinates. To be specific and to keep things consistent with the rest of my project, I'm using hg38 on the website genome.ucsc.edu. So far, I've been typing one gene symbol at a time into the gene browser and clicking the correct match. This sends me to the gene's page, which shows me the start and end genome coordinates for the gene's transcript and coding region. That's what I'm looking for.
This is a tedious process when I have hundreds of genes to go through, so I was wondering if I can input a list and get out a list instead. I know about the page https://genome.ucsc.edu/cgi-bin/hgTables. As an example, suppose I click "paste list" and just type the following two lines:
CASZ1
MTHFR

Here's a screenshot of the output (click to see full sized image):

I think this page is generally intended to go the opposite direction (genomic coordinates -> genes)?
Edit: Ok, here's an example of what I want. I'd like to just enter the gene symbols and get the transcript and coding regions shown.
CASZ1, Transcript chr1:10,636,604-10,796,646, Coding Region chr1:10,638,942-10,693,889
MTHFR, Transcript chr1:11,785,723-11,805,964, Coding Region chr1:11,790,680-11,803,116
ST7L, Transcript chr1:112,525,883-112,618,928, Coding Region chr1:112,526,013-112,598,043

I obtained these by searching through genome.ucsc.edu one by one manually.

Does it have to be from the UCSC browser? It isn't the best tool for this sort of thing. Are you open to other tools?

@terdon I prefer UCSC browser because I used it for gene mapping earlier in this project. Redoing everything after that step is too costly at this point. So if I switch to another website, I'll probably have to accept some inconsistency. For instance, I tried looking up a few genomic coordinates on https://www.ncbi.nlm.nih.gov/datasets/tables/genes/. Coordinates are roughly similar to UCSC's transcript coordinates, but a little bit different.
Unless there's a tool that can recover exactly what UCSC would have given me?

Yes, you can get the data directly from the UCSC database. I just wouldn't use the browser for this. Can you give us some input data that we can use to test any solution? Just [edit] your question and show us a few lines of what you have (different IDs would require different approaches). Also, please clarify if you want the location of the gene (depending on your definition, this can include promoters and other regulatory sequences) or just the gene's various transcripts.

Does this answer your question? Create BED file from list with gene symbols

@terdon Ok, I have updated the question. Just taking a quick look at user324810 link though, it seems to assume familiarity with something that I haven't seen before; so I'll have to look into how to use that. Thanks

@user324810 Regarding your answer in that question link you posted, yes, that's the webpage I tried. As I wrote in my question though, it returns a long output page when I input just 2 symbols, which seems not right. Even stranger is why it provides the option of specifying gene position, when the goal is to find the gene position!
I'll add some more clarification to my question.

The current release of Gencode for hg38 is v38. This is one way you could get a list of Gencode gene annotations:
$ wget -qO- ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_38/gencode.v38.annotation.gff3.gz \
    | gunzip --stdout - \
    | awk '$3 == "gene"' - \
    | convert2bed -i gff - \
    > genes.bed

To get convert2bed, you could install a toolkit called BEDOPS, which is available via Github or package managers. Refer to the installation instructions for more detail.
Once you have genes.bed, you don't need to make it again, unless Gencode publishes a new set of annotations.
You could create a text file called genes_of_interest.txt that contains your HGNC names of interest, e.g.:
$ cat > genes_of_interest.txt
CASZ1
MTHFR
ST7L

Press Ctrl + D at the end to close the file.
To search genes.bed for entries where there are matches to your genes of interest, you can use grep:
$ grep -Fwf genes_of_interest.txt genes.bed
chr1    10636603    10796650    ENSG00000130940.14  .   -   HAVANA  gene    .   ID=ENSG00000130940.14;gene_id=ENSG00000130940.14;gene_type=protein_coding;gene_name=CASZ1;level=1;havana_gene=OTTHUMG00000002035.5
chr1    11785722    11806920    ENSG00000177000.11  .   -   HAVANA  gene    .   ID=ENSG00000177000.11;gene_id=ENSG00000177000.11;gene_type=protein_coding;gene_name=MTHFR;level=2;havana_gene=OTTHUMG00000002277.5
chr1    112523517   112620825   ENSG00000007341.18  .   -   HAVANA  gene    .   ID=ENSG00000007341.18;gene_id=ENSG00000007341.18;gene_type=protein_coding;gene_name=ST7L;level=2;havana_gene=OTTHUMG00000011753.7

The above is just sample output. Your result may differ.


text:
Creating Experimental Components on SBOL
I am developing software to store experiments or studies using SBOL3. In the regular synthetic biology workflow we will find the experimental components Strain and Media. A safe way to implement them is using tyto.SBO.functional_entity as type for both.
Other attempts to implement media (https://github.com/Bioprotocols/paml/blob/main/examples/growth_curve.py)looks like this:
SC_media = sbol3.Component('SC_Media', 'TBD', name='Synthetic Complete Media')

Therefore there is a need for standardise these types on SBOL3 to provide more accurate representations. The Experimental Factor Ontology (EFO) have ontology terms for Strain http://www.ebi.ac.uk/efo/EFO_0005135 and Media http://www.ebi.ac.uk/efo/EFO_0000579 that can be used as type. Then draft SBOL components looks as follow:
strain = sbol3.Component(fj_strain.name(), http://www.ebi.ac.uk/efo/EFO_0005135)
strain.name =fj_strain.name
strain.description = fj_strain.description()

media = sbol3.Component(fj_media.name(), http://www.ebi.ac.uk/efo/EFO_0000579)
media.name =fj_media.name
media.description = fj_media.description()

The problem with this is the lack of readability on the code given the obscure ontology URIs. This problem has been solved for other ontologies using tyto, but EFO is not supported out-of-the-box and the configuration of new ontologies is not working. On tyto repository there is an example:
from tyto import EBIOntologyLookupService, Ontology
KISAO = Ontology(uri='http://www.biomodels.net/kisao/KISAO_FULL#', endpoints=[EBIOntologyLookupService])
KISAO.Gillespie_direct_algorithm
'http://www.biomodels.net/kisao/KISAO#KISAO_0000029'

But I copied and runed that code but it raise this error:
[ERROR] tyto.py 70: Ontology http://www.biomodels.net/kisao/KISAO_FULL# is not available at EBI Ontology Lookup Service
What I would like to do is something like this:
from tyto import EBIOntologyLookupService, Ontology
efo = Ontology(uri='http://www.ebi.ac.uk/efo', endpoints=[EBIOntologyLookupService])
efo.media

But this raises the same error.
How can I add EFO to tyto? What is the URI that should be provided? Should I raise an issue in tyto to ask for EFO out-of-the-box support?
Do you agree with the ontology terms for Strain and Media? Any feedback?

Section 7.9 of the SBOL3 specification already provides recommended terms for experimental elements.  Among these, Section 7.9.2 specifically addresses strain and media, reading in pat:

The roles of reagents may vary in context: for example, arabinose may serve as an inducer or as a media carbon source. As such, contextual role SHOULD be indicated by an NCI Thesaurus (NCIT) term in a role property of the Feature. Examples include:

NCIT:C64356: PositiveControl
NCIT:C48694: Cell
NCIT:C85504: Media
NCIT:C14419: Strain
NCIT:C120268: Inducer

NCIT is one of the ontologies that is built-in for tyto:
tyto.NCIT.Strain
Out[4]: 'https://identifiers.org/ncit:C95501'
tyto.NCIT.Media
Out[5]: 'https://identifiers.org/ncit:C19442'

One thing that makes noise to me is that these terms are exemplified as roles not as types. These terms make sense as types to me.  Should I use another type like SBO_FUNCTIONAL_ENTITY and keep the NCIT terms as roles?

Type is used to indicate the nature of the material. Most of these are complex entities that cannot be categorized as being a single class of molecule, and thus SBO_FUNCTIONAL_ENTITY is the appropriate type.


text:
How to compare sequences of genes obtained after whole-genome sequencing to reference genome?
My idea is to align whole-genome sequencing data (as fastq files after, 30× coverage, gDNA) to the reference mouse genome (NCBI), extract the immunoglobulin locus and compare it to the reference. I think the question is too broad but maybe a hint how to start ? (I'm new to the field.)
Task is: to check for difference between reference genome and a sample for 100-200 genes.
I checked these: question (but it's closed), answer (gave a hint about whole reference mapping not only the locus of interest), question (not answered), question (useful I guess but not sure).
Additional info:
Instrument Platform: ILLUMINA
Instrument Model: HiSeq X Ten
Library Layout: PAIRED
Library Strategy: WGS
Library Source: GENOMIC
Typical sequence length distribution from FastQC: 150nt.
I have Win10 with WSL 1 installed.

Please [edit] your question and give us more details. What species are you working with? What technique? Is this DNAseq or RNAseq or something else? Do you know that the locus is or is that what you want to find out?

@terdon: Thanks for the comment. I've edited. But I didn't understand your question "What technique?". What do you mean?

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

What have you tried?

I've managed to get BAM and BAI files with https://usegalaxy.org. Now, I'm trying to visualize with IGV...

Surprisingly there aren't actually many helpful one-liners to map Illumina WGS reads on the internet - everyone just sort of assumes that people can figure it out from the manual. Many people end up doing it in a not-so-ideal way, in that they stop at the .sam files that bwa mem makes that take up a lot of disk space.
I like this incantation:
# index the reference
bwa index reference_genome.fasta

# map your reads to the reference (bwa mem)
# compress that output to a bam file (samtools view -hb)
# sort the bam file (samtools sort)

bwa mem -t <number_of_threads> reference_genome.fasta your_sample.R1.fastq.gz your_sample.R2.fastq.gz | \
samtools view -hb - | \
samtools sort - > your_reads_to_reference.sorted.bam

# index the bamfile. makes a .bai file

samtools index your_reads_to_reference.sorted.bam

With this sorted .bam, the .bai, and the reference genome you can now use IGV to directly look at the reads mapped to your genome.
Also with the bam, you can generate a VCF file using freebayes or another tool like samtools mpileup. This will give you a file of how your sample is different from the reference genome. You can use bedtools intersect to get the variants only from your gene/regions of choice.

conchoecia, thank you for your answer and suggestions. I did Bowtie2 at https://usegalaxy.org. Liam McIntyre and you suggest bwa mem, minimap2, or ngmlr. My reads are ~150nt. What tool should I choose? I guess for 150nt bwa mem is better. Currently with .bam and .bai files I'm trying to visualize data in IGV, but yes, genes of interest are preferred to show. I hope to completely switch to command line

bwa mem is the "standard" for Illumina WGS reads. minimap2 is the "standard" for long WGS reads, like PacBio or ONT.

Then I'll do bwa mem. Thank you

according to timpo0's comment https://github.com/jts/methylation-analysis/issues/13#issuecomment-370628831 with reference to  https://github.com/lh3/bwa ("For all the algorithms, BWA first needs to construct the FM-index for the reference genome (the index command)."): bwa index reference.fasta precede the rest lines

.bam file's size is 0 bytes...

@abc I also got 0 bytes also.  here is a possible helper bwa index ref.fa; bwa mem ref.fa out1.fq out2.fq |  samtools sort - -o out.sorted.bam. I am using samtools 1.14 and bwa-0.7.17-r1188

you can also try this command from the samtools docs "fastq to bam", http://www.htslib.org/workflow/#fastq_to_bam just swap minimap2 with bwa mem
e.g. bwa mem ref.fa out1.fq out2.fq | samtools fixmate -u -m - - | samtools sort -u -@2 - | samtools markdup -@8 --reference ref.fa - final.cram

@ColinD: thanks for suggestion, 9h run (1st code) resulted in 12 out.sorted.bam.tmp.00XX.bam files of ~260Mb, the run is still going. A file is generated ~ every 30mins. In shell I see [M::process] read 70190 sequences (10000248 bp)... [M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (8, 26835, 0, 4) [] skip orientation FF as there are not enough pairs [] analyzing insert size distribution for orientation FR... [] (25, 50, 75) percentile: (172, 245, 331) [] low and high boundaries for computing mean and std.dev: (1, 649) [] mean and std.dev: (257.71, 118.14) <...>

<...> [] low and high boundaries for proper pairs: (1, 808) [] skip orientation RF as there are not enough pairs [] skip orientation RR as there are not enough pairs [M::mem_process_seqs] Processed 70060 reads in 64.900 CPU sec, 64.655 real sec

that all sounds expected :)

This is very easy. First map all reads to the reference sequence for the organism. If they are short reads use BWA mem. If long reads use minimap2 or ngmlr. Once mapped view the locus of interest in igv (GUI). If you need to extract the reads at the locus pass the coordinates to samtools view (cli).

MvIntyre: thank you for the answer, may be you can give advice any source on the topic as well ? I found BWA-MEM: https://github.com/lh3/bwa and https://www.ridom.de/u/BWA_Mapper.html, IGV: http://software.broadinstitute.org/software/igv/home, Samtools: http://www.htslib.org/

They are all git projects


text:
How to identify Father-Daughter Relationship from Phylogenetic analysis?
I have performed a phylogenetic analysis using CLUSTALW and obtained a phylogenetic tree. I know how to read a phylogenetic tree and how to relate the sequence similarity between two sequences by looking at the tree, but what I am concerned about is how one can decipher the father-daughter relationship by looking at the tree? Also, if I want to isolate the sequence of the father corresponding to a daughter sequence how can I do that?
I have attached the corresponding image. In the image, I know that sequences "QXN18196" and "QXN18436" are closely related. But will it be safe to say that sequence "QXN18520" is the father sequence for both the sequences? Or if I want to establish a father-daughter relationship for the whole tree, what will be the relation between the sequences?

Better terminology in this context is to use Parent-Children/child relationships, rather than Father-Daughter

Thank you so much for your answer. I am still not sure how can one establish that relation from phylogeny. Or to put it another way, what is the right approach to determine a parent-child relationship between a set of sequences, say 6 sequences.

Unless you have wildly differing mutation rates on different branches, you can't say that QXN18520 has the same sequence as the most recent common ancestor of QXN18196 and QXN18436, because QXN18520 was collected today, not back whenever QXN18196 and QXN18436 shared a common ancestor. All that you can say from this phylogeny is that there probably existed a common ancestor of QXN18196 and QXN18436, marked here by 5e-9[86.8], and that this common ancestor and QXN18520 shared a common ancestor even further back, marked here by 5e-9[0].
(General observation: Note also that for a number of reasons, phylogenies often don't reflect reality -- they are simply a maximum likelihood estimate of the actual ancestral pedigree, often under some very unrealistic assumptions, and/or according to an incomplete or incorrect view of reality -- e.g. the samples that you have collected may have a lot of allele matches due to identity by state (IBS) rather than identity by descent (IBD). In other words, let phylogenies inform you, but don't assume they are the absolute truth.)

Thank you so much for your answer. I was wondering then how can one establish a father-daughter relation from the phylogenetic tree? If phylogeny is an estimate of ancestral pedigree, then how can one find out if a set of sequences really do have a father-daughter relationship?

@satyamsangeet if you suspect there is a literal father and a literal daughter in the data somewhere, and you want to detect this, then you first need to know which chromosome the alleles came from. For the autosomal diploid chromosomes (1-22), inheritance is via random assortment, and you would expect very close to 50% of the alleles will match between father and daughter, but only if all markers are widely spaced apart (so that you have low linkage disequilibrium). If all or even some alleles are in high linkage disequilibrium, then allele sharing can deviate from 50%. (ctd...)

However what I just said about 50% sharing between father and daughter is only really true if you have diploid data for both individuals (meaning you have the alleles for a given locus from both chromosomes). If you have alleles from only one chromosome for the father and the daughter, then you would expect only 25% sequence identity (sequence match), because you're missing half the data. Also, note that you can't tell who the father is and who the daughter is simply by looking at the percentage of sequence identity.

@satyamsangeet Chromosome X is inherited differently to chromosomes 1-22 between father and child. You would expect a 100% match between the father's alleles and one of the daughter's alleles, if you have diploid data for the daughter. The father has only one copy of the X chromosome, so you would only ever have haploid data for the father, and you would be able to tell he's the male based on this, if in general your data is diploid. If you have only haploid data for all samples (if you have alleles from only one of the daughter's X chromosomes), then you would expect sequence identity of 50%.

Then the Y chromosome is of course not even possessed by the daughter, but will be present for the father if it has been sampled (this is another way to tell the father and daughter apart). And the mitochondrial chromosome will be present for all individuals, if it has been sampled, but it will not match between a father and a daughter in general -- it should only match between a daughter and her mother (or a son and his mother).


text:
How to train annotations tools
I would like to train Augustus, SNAP and GlimmerHMM. I found protein sequences in GenBank and orthodb.org. Furthermore, I found HMM files on busco-data.ezlab.org.
$ wget -c https://busco-data.ezlab.org/v5/data/lineages/viridiplantae_odb10.2020-09-10.tar.gz
$ wget -c https://v100.orthodb.org/download/odb10_plants_fasta.tar.gz

Are there any instructions on how to train those annotations tools?
Thank you in advance


text:
Organelle genome genes
How can I get all proteins coded by mitochondrial genome or plastid genome from UniProt. I see UniProt has multiple filtering and search parameters but I could not distinguish between proteins/genes coded by nuclear DNA and mitochondrial or plastid genomes.
I specifically want UniProt (rather than e.g. Organelle Genome Resources) because usually annotation of proteins uses UniProt ids.
Any suggestions?

You can use the UniProt Advanced search for this.
If you look at this help page https://www.uniprot.org/help/advanced_search, in particular the first example screenshot, you can type "encoded" where the screenshot has "structure". Then you can choose your organelle. For Plastid, for example, this will translate into the query syntax organelle:plastid :
https://www.uniprot.org/uniprot/?query=organelle%3APlastid&sort=score
Several other organelles are available, of course, e.g. Mitochondrion :
https://www.uniprot.org/uniprot/?query=organelle%3AMitochondrion&sort=score
Please don't hesitate to contact the UniProt helpdesk. We can unfortunately not guarantee to monitor all available forums as closely as we would like to, and can answer much more promptly when contacted directly.


text:
How to de novo hybrid assemble with Pacbio CCS and Illumina PE reads
I would like to perform de novo genome assembly on a diploid microalgal strain.
I have two datasets:

PacBio CCS/HiFi reads, low coverage.
Illumina PE 2x150 (standard shotgun)

Does anybody have any software recommendations for how best to use both in a hybrid de novo assembly?
Most options out there are set up for Pacbio CLR reads, not Pacbio CCS/HiFi, for using in hybrid assembly. Perhaps I should assemble short read contigs first then combine with wengan, or go back to the PacBio subreads instead of CCS to use in a hybrid method?

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

If you have reasonable coverage with HiFi (30X per haplotype), then there usually isn't any reason to also use Illumina. My understanding is that Illumina can make the assembly worse if you try to polish with it, because it can erroneously induce haplotype switching errors.
I would use HiCanu or HiFiAsm, depending on whether you expect significant heterozygosity, and also assuming that you have non-trivial HiFi coverage.
If your HiFi coverage is very low (<15X per haplotype, maybe), you might consider going back to subreads, assembling those, and polishing with Pilon and your Illumina reads, or one of the workflows that you have found. But otherwise I don't think the Illumina will be very useful. If you had Hi-C reads you could use them directly in HiFiAsm to assist phasing and scaffolding, but I'm assuming that your Illumina is shotgun.
I know it's a PB promo material, but you see that PB promotes HiFi as not needing polishing, as HiFi has accuracy >99%, so shotgun Illumina is not going to do much. There is some argument for still polishing with your original long reads.
I can recommend any of the linked papers for ideas or further explanation.

Unfortunately, my HiFi coverage isn't that great, (unexpectedly large genome). However, my shotgun short read dataset is large and good quality, hence the question.

@bishopia Hmmmm I see what you mean. In that case, subread assembly+Illumina polishing might be your best bet. It would still be helpful to know more specifics (e.g. coverage, read length profile of each dataset: CCS, subread, ILMN). I am pretty skeptical of the usefulness of Illumina contigs as the basis of assemblies of large/complex/heterozygous genomes, compared to even low-coverage PB data.


text:
How can I get a list of neighboring Hydrogen atoms of an alpha-carbon?
Suppose, I want to find the neighboring Hydrogen atoms of a specific C-alpha atom.
The following was my effort. However, it produces a bunch of empty lists.
from typing import List
from Bio.PDB import Chain, Atom, NeighborSearch, PDBParser

def get_hydrogen_atoms(c_alpha_atom):
    ns = NeighborSearch([c_alpha_atom])
    return ns.search_all(190.0, "A")

p = PDBParser()
structure = p.get_structure('16pk', '16pk.pdb')

c_alpha_list = [atom for atom in structure.get_atoms()]

for c_alpha in c_alpha_list:
    if c_alpha.name == "CA":
        hydrogen_list = get_hydrogen_atoms(c_alpha)
        print(hydrogen_list)

How can I correct this?

hydrogen_list = get_hydrogen_atoms(c_alpha) should be hydrogen_list.append(get_hydrogen_atoms(c_alpha)).
Also your structure needs to be protonated —most crystal structures are not. radius is in Å, not pm. Cα atoms are no hydrogen bond donors/acceptors. Why no search for nearby N atomnames to the atomname O, which are donor - acceptor pairs...

You posted similar question here: https://stackoverflow.com/questions/69676943/how-can-i-get-a-list-of-neighboring-hydrogen-atoms-of-an-alpha-carbon/69815509#69815509
ok got some code :
from Bio.PDB import Chain, Atom, NeighborSearch, PDBParser, Selection

def get_hydrogen_atoms(c_alpha_atom: Atom.Atom, po):
    ns = NeighborSearch(po)
    return [h for h in ns.search( c_alpha_atom.get_coord(), 10.0, 'A') if h.element == 'H']

p = PDBParser()
# structure = p.get_structure('16pk', '16pk.pdb')
structure = p.get_structure('16pk', 'pH8_5_1bq9_short.pdb')

c_alpha_list = [atom for atom in structure.get_atoms() if atom.name == "CA"]

print('c_alpha_list[0]  : ', c_alpha_list[0], c_alpha_list[0].id , c_alpha_list[0].get_serial_number(),'\n')

po = Selection.unfold_entities(structure, 'A')

hydrogen_list = get_hydrogen_atoms(c_alpha_list[0], po)

print('hydrogen_list : \n',hydrogen_list)

print('\nSelection (po) : \n', po)

had to use a different PDB file because I needed the H atoms,
pH8_5_1bq9_short.pdb :
HETATM    1  N   FME A   1      23.447  -5.765   5.119  1.00 19.07           N
HETATM    2  CN  FME A   1      23.308  -6.776   4.291  1.00 20.46           C
HETATM    3  O1  FME A   1      22.504  -6.713   3.357  1.00 24.75           O
HETATM    4  CA  FME A   1      22.488  -4.745   5.410  1.00 18.15           C
HETATM    5  CB  FME A   1      22.433  -4.342   6.827  1.00 21.15           C
HETATM    6  CG  FME A   1      22.198  -5.618   7.751  1.00 22.84           C
HETATM    7  SD  FME A   1      20.445  -6.032   7.772  1.00 23.40           S
HETATM    8  CE  FME A   1      20.476  -7.587   8.656  1.00 28.76           C
HETATM    9  C   FME A   1      22.699  -3.530   4.487  1.00 17.12           C
HETATM   10  O   FME A   1      23.804  -3.248   4.077  1.00 20.25           O
HETATM    0 3HE  FME A   1      19.574  -7.938   8.726  1.00 28.76           H   new
HETATM    0 2HG  FME A   1      22.713  -6.369   7.418  1.00 22.84           H   new
HETATM    0 2HE  FME A   1      21.032  -8.221   8.177  1.00 28.76           H   new
HETATM    0 2HB  FME A   1      23.260  -3.901   7.077  1.00 21.15           H   new
HETATM    0 1HG  FME A   1      22.508  -5.438   8.652  1.00 22.84           H   new
HETATM    0 1HE  FME A   1      20.838  -7.447   9.545  1.00 28.76           H   new
HETATM    0 1HB  FME A   1      21.717  -3.701   6.961  1.00 21.15           H   new
HETATM    0  HN  FME A   1      24.199  -5.711   5.532  1.00 19.07           H   new
HETATM    0  HCN FME A   1      23.841  -7.577   4.415  1.00 20.46           H   new
HETATM    0  HA  FME A   1      21.619  -5.138   5.231  1.00 18.15           H   new
ATOM     11  N   ALA A   2      21.633  -2.939   4.043  1.00 15.50           N
ATOM     12  CA  ALA A   2      21.644  -1.884   3.062  1.00 15.08           C
ATOM     13  C   ALA A   2      20.899  -0.661   3.519  1.00 13.72           C
ATOM     14  O   ALA A   2      20.278  -0.702   4.584  1.00 13.49           O
ATOM     15  CB  ALA A   2      21.083  -2.339   1.710  1.00 17.86           C
ATOM      0  H   ALA A   2      20.842  -3.144   4.312  1.00 15.50           H   new
ATOM      0  HA  ALA A   2      22.579  -1.652   2.951  1.00 15.08           H   new
ATOM      0  HB1 ALA A   2      21.110  -1.599   1.083  1.00 17.86           H   new
ATOM      0  HB2 ALA A   2      21.618  -3.072   1.368  1.00 17.86           H   new
ATOM      0  HB3 ALA A   2      20.166  -2.634   1.823  1.00 17.86           H   new
ATOM     16  N   LYS A   3      21.005   0.412   2.774  1.00 12.48           N
ATOM     17  CA  LYS A   3      20.341   1.658   3.029  1.00 12.61           C
ATOM     18  C   LYS A   3      19.415   1.991   1.879  1.00 11.53           C
ATOM     19  O   LYS A   3      19.784   1.753   0.701  1.00 11.66           O
ATOM     20  CB  LYS A   3      21.338   2.806   3.273  1.00 13.38           C
ATOM     21  CG  LYS A   3      22.062   2.721   4.593  1.00 16.70           C
ATOM     22  CD  LYS A   3      23.062   3.847   4.814  1.00 18.05           C
ATOM     23  CE  LYS A   3      23.806   3.576   6.136  1.00 20.42           C
ATOM     24  NZ  LYS A   3      24.878   4.503   6.403  1.00 22.91           N
ATOM      0  H   LYS A   3      21.494   0.433   2.067  1.00 12.48           H   new
ATOM      0  HA  LYS A   3      19.821   1.557   3.842  1.00 12.61           H   new
ATOM      0  HB2 LYS A   3      21.992   2.812   2.556  1.00 13.38           H   new
ATOM      0  HB3 LYS A   3      20.862   3.650   3.230  1.00 13.38           H   new
ATOM      0  HG2 LYS A   3      21.411   2.732   5.312  1.00 16.70           H   new
ATOM      0  HG3 LYS A   3      22.527   1.871   4.644  1.00 16.70           H   new
ATOM      0  HD2 LYS A   3      23.689   3.892   4.076  1.00 18.05           H   new
ATOM      0  HD3 LYS A   3      22.606   4.702   4.851  1.00 18.05           H   new
ATOM      0  HE2 LYS A   3      23.170   3.611   6.868  1.00 20.42           H   new
ATOM      0  HE3 LYS A   3      24.166   2.676   6.116  1.00 20.42           H   new
ATOM      0  HZ1 LYS A   3      25.266   4.294   7.176  1.00 22.91           H   new
ATOM      0  HZ2 LYS A   3      25.479   4.458   5.748  1.00 22.91           H   new
ATOM      0  HZ3 LYS A   3      24.551   5.330   6.449  1.00 22.91           H   new
ATOM     25  N   TRP A   4      18.274   2.598   2.165  1.00 10.75           N
ATOM     26  CA  TRP A   4      17.283   3.019   1.181  1.00 10.69           C
ATOM     27  C   TRP A   4      16.875   4.472   1.502  1.00 11.37           C
ATOM     28  O   TRP A   4      16.812   4.872   2.651  1.00 17.19           O
ATOM     29  CB  TRP A   4      16.012   2.143   1.352  1.00 11.20           C
ATOM     30  CG  TRP A   4      16.164   0.734   0.997  1.00 11.80           C
ATOM     31  CD1 TRP A   4      16.913  -0.218   1.669  1.00 12.50           C
ATOM     32  CD2 TRP A   4      15.535   0.036  -0.079  1.00 11.53           C
ATOM     33  NE1 TRP A   4      16.806  -1.435   1.050  1.00 12.76           N
ATOM     34  CE2 TRP A   4      15.959  -1.302  -0.013  1.00 12.61           C
ATOM     35  CE3 TRP A   4      14.642   0.421  -1.079  1.00 11.88           C
ATOM     36  CZ2 TRP A   4      15.531  -2.258  -0.929  1.00 13.18           C
ATOM     37  CZ3 TRP A   4      14.217  -0.527  -1.989  1.00 12.54           C
ATOM     38  CH2 TRP A   4      14.668  -1.858  -1.894  1.00 14.26           C

output :
c_alpha_list[0]  :  <Atom CA> CA 4

hydrogen_list :
 [<Atom HB1>, <Atom HB3>, <Atom H>, <Atom 2HE>, <Atom 3HE>, <Atom 1HE>, <Atom HB3>, <Atom HA>, <Atom HG2>, <Atom HB2>, <Atom HCN>, <Atom HA>, <Atom HN>, <Atom 1HB>, <Atom 2HB>, <Atom 2HG>, <Atom 1HG>, <Atom H>, <Atom HB2>, <Atom HA>, <Atom HD2>, <Atom HG3>, <Atom HD3>, <Atom HZ2>, <Atom HE3>, <Atom HE2>, <Atom HZ1>]

Selection (po) :
 [<Atom N>, <Atom CN>, <Atom O1>, <Atom CA>, <Atom CB>, <Atom CG>, <Atom SD>, <Atom CE>, <Atom C>, <Atom O>, <Atom 3HE>, <Atom 2HG>, <Atom 2HE>, <Atom 2HB>, <Atom 1HG>, <Atom 1HE>, <Atom 1HB>, <Atom HN>, <Atom HCN>, <Atom HA>, <Atom N>, <Atom CA>, <Atom C>, <Atom O>, <Atom CB>, <Atom H>, <Atom HA>, <Atom HB1>, <Atom HB2>, <Atom HB3>, <Atom N>, <Atom CA>, <Atom C>, <Atom O>, <Atom CB>, <Atom CG>, <Atom CD>, <Atom CE>, <Atom NZ>, <Atom H>, <Atom HA>, <Atom HB2>, <Atom HB3>, <Atom HG2>, <Atom HG3>, <Atom HD2>, <Atom HD3>, <Atom HE2>, <Atom HE3>, <Atom HZ1>, <Atom HZ2>, <Atom HZ3>, <Atom N>, <Atom CA>, <Atom C>, <Atom O>, <Atom CB>, <Atom CG>, <Atom CD1>, <Atom CD2>, <Atom NE1>, <Atom CE2>, <Atom CE3>, <Atom CZ2>, <Atom CZ3>, <Atom CH2>]

the selection list contains both ATOM and HETATM
have a look at : https://www.biostars.org/p/269579/ How to use biopython neighborsearch?
and Bio.PDB.NeighborSearch module :https://biopython.org/docs/1.75/api/Bio.PDB.NeighborSearch.html ;
Bio.PDB.Atom module :https://biopython.org/docs/1.75/api/Bio.PDB.Atom.html


text:
Error in "apply(grps, 1, function(x) any(is.na(x))) : dim(X) must have a positive length" at running dht2 function in R
I am running the script from de Distance examples (http://examples.distancesampling.org/Distance-cues/cuecounts-distill.html) to estimate density using function dht2 of Distance library but I have an error.
The script is:
library (Distance)
data(wren_cuecount)
conversion.factor <- convert_units("meter", NULL, "hectare")
wren_cuecount$Effort <- wren_cuecount$Search.time
wrensong.hr <- ds(wren_cuecount, transect="point", key="hr", adjustment=NULL,
                  truncation=100)
cuerate <- unique(wren_cuecount[ , c("Cue.rate","Cue.rate.SE")])
names(cuerate) <- c("rate", "SE")
(mult <- list(creation=cuerate))
wren.estimate <- dht2(wrensong.hr, flatfile=wren_cuecount, strat_formula=~1,
                      multipliers=mult, convert_units=conversion.factor)
print(wren.estimate, report="density")""

When I run the dht2 function I get:

Error in apply(grps, 1, function(x) any(is.na(x))) :    dim(X) must
have a positive length

How can I solve this?

Please format your code so that it can be understood.

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

I would suggest inspecting some of the intermediate objects such as wrengsong.hr, wren_cuecount, mult, etc. to ensure that they match the expected inputs of those functions.


text:
How can I programmatically add Hydrogen to a PDB structure using BioPython?
I need to programmatically add Hydrogen atoms to a PDB file (say, 16PK.pdb) using BioPython.
What algorithm should I use for this?

Not actually and answer as you cannot do it with biopython.PDB. But there are may other tools.
Several are online, for example: https://zhanggroup.org/HAAD/ by the folk who wrote I-Tasser.
Nearly every MD pipeline has their own protonation algorithm. The only hydrogens that are an issue are histidine protonation —HIE vs. HID— and protonation states of catalytic residues (CYZ vs. CYS etc).
The PyMOL python module is really good and can do a lot. Only issue is that the easiest way to install is via conda.
conda install -c conda-forge -c schrodinger pymol-bundle

with which one can:
import pymol2
pymol2.PyMOL() as pymol:
    pymol.cmd.load('.pdb', 'myprotein') # file
    # pymol.cmd.fetch('', 'myprotein') # PDB code
    pymol.cmd.h_add()
    pymol.cmd.save('.pdb')

One can find neighbours to a given atom too (see this Q&A), select by secondary structure, define secondary structure (pymol.cmd.dss()) etc. A caveat is that, bar for some commands like pymol.cmd.get_model(selection).atom, one does not interact directly with the atoms and everything is done via selections.
However, as this Q follows on from this unanswerable one and in turn from several others about SS. I am guessing it is to determine the hydrogen bond distances between different sheets. These have a backbone amide nitrogen (name N) donating a hydrogen (H, not HA etc.) to an accepting backbone amide oxygen of another sheet (O). Unless you want an exact distance, you could simply use the distance between these heavy atoms, which would be slightly larger than the H-N bond plus the H•O hydrogen bond and the oxygen may not be coplanar with the planar amide.

thanks a lot...

Would it be a realistic goal to implement HAAD protonation algorithm with BioPython?

Doable totally, worth doing something that is implemented in many places elsewhere on a module that's a swiss-army knife, but limited in advanced functionality, definitely no...

The incorrect syntax pymol2.PyMOL() as pymol: should throw an error.

Not on python 3. In this case it's not the module pymol but an instance of multisession pymol2.PyMOL which behaves like the module but is way better

conda install -c conda-forge openbabel
from openbabel import pybel as pb
babel_mol = pb.readfile(format="pdb", filename=path_to_pdb).__next__()
babel_mol.addh()
babel_mol.write(format='pdb', filename=path_to_babel_pdbh, overwrite=True)


text:
Metagene analysis with deepTools; 5UTR-CDS-3UTR and Exon-Intron-Exon
I have some PCR-CLIP data. PAR-CLIP identifies targets for RNA binding proteins. I am interested to build metagene profile so that transcript wide binding pattern can be compared. What I am mostly interested to plot:

5UTR-CDS-3UTR
Exon-Intron-Exon

With deepTools I can make a combined profile where metagene profile for whole transcript length is found.
For 5UTR-CDS-3UTR, what I did:

create bigwig files (from bam files) for signal

used GTF as region of interest

The code goes as follows:
computeMatrix scale-regions
--metagene
-m 1000
-S *bw
-R *gtf
-a 1000
-b 1000
--transcriptID transcript
--transcript_id_designator transcript_id
-o outMat.mtx

plotProfile -m outMat.mtx -out outFig.png --perGroup

Generated image, given below, where my understanding is

-1.0:TSS is 5' UTR,
TSS:TES is CDS+3'UTR

If that is so, I am missing 3' UTR region.
So what I did, running computeMatrix reference-point by centering at TSS, center and TES with +/- 1000bp. The plot is as follows:

With some operations I converted it to:

Beside this kind of manual fix, I believe there are some systematic way to render such metagene profile where 5UTR-CDS-3UTR would be side by side for comparison.
1. For plotting such 5UTR-CDS-3UTR metagene profile systematically, how should I proceed with deepTools
In addition to that if I am interested to plot such metagene profile for 5Exon-Intron-3Exon, like the example given below.

2. For generating Exon-Intron-Exon metagene profile how could I proceed with deepTools?


text:
Do I have to download all nt.xx files to perform a search in blastn 2.12.0
I downloaded and extract blast nt.24 database, but when I run :
blastn -task blastn-short -query consulta.fa -db nt.24 -out misalida.txt
it says: mdb_env_open: No such file or directory
and if a try with nt instead of nt.24 it says:
BLAST Database error: Could not find volume or alias file (nt.00) referenced in alias file (/root/my-documents/blastdb/nt/nt
So my question is :
Do I have to download all nt.xx files?
And : Will it show the organisms in the output result?

In general, yes, you do need to download and extract all (nt.00 - nt.46 as of now) of the nt portions in order to do a comprehensive BLAST against all of nt. The database has been split up for convenience (to make smaller individual files), so all are necessary to be comprehensive.
Once you download all the pieces, you’ll want to use nt as the database parameter to blastn.
blastn -task blastn-short -query consulta.fa -db nt -out misalida.txt
As for the second question - this depends on how you run the BLAST. Each accession is associated with a taxid, which will map to its taxonomy. This can be mapped back to the organism, but you may need to do some interpretation in order to convert the BLAST hits to the proper organism - this depends highly on what kind of sequences you are aligning.

Thanks, I downloaded nt.00 (by http firefox because wget via ftp always returned a corrupted file when decompressing it) and it works fine locally.
The command is:
blastn -task blastn-short -query consulta.fa -db nt.00 -out misalida.txt
I also added a first line >Seq1 to the query file consulta.fa (if not didn't work):>Seq1
GTGCC...
And the results in misalida.txt first print the taxid with the organism name and chromosome. At the bottom it printed the positions details.
What I'm doing is a genome editor, so when you enter some nucleotids(like 14 to 20 or more) it will search with blastn possible next nucleotids paths from many organisms with  blastdbcmd.


text:
In Seurat, how can I export cluster IDs to csv files
I am doing scRNAseq analysis with Seurat.
I clustered the cells using the FindClusters() function.
What I want to do is to export information about which cells belong to which clusters to a CSV file.
In a Seurat object, we can show the cluster IDs by using Idents(・), but I have no idea how to export this to CSV files.
I would be grateful if you could show this by using the PMBC data (https://satijalab.org/seurat/articles/pbmc3k_tutorial.html) as an example.

The clusters information, together with additional information on each cells (%mito, sample name, group, orig.ident,...) are stored inside seuratObject@meta.data.
write.csv(seuratObject@meta.data,"./seurat_metadata.csv") will dump everything to a file.
The method @osmoc suggested, is using built in seurat functions to interact with seuratObject@meta.data.
seuratObject$graph_name_resolution points to seuratObject@meta.data$graph_name_resolution.
I personally prefer to use the metadata table, as I have more control on how I interact with the data.

Thank you very much !

Using the meta.data slot is the best answer. The behavior of various shortcut methods can change between different versions of Seurat and even based on how you interact with your Seurat object. More generally, it is a bad idea to take columns from different sources when constructing a data frame.

I've recently done this with my own Seurat data; here's the process:
# load tidyverse packages
library(tidyverse)
# create cluster ID table
tibble(cellID = colnames(pbmc), clusterID = Idents(pbmc)) %>%
  # write out to a file, with today's date
  write_csv(file = sprintf("cluster_mappings_%s.csv", Sys.Date()))

Thank you very much !!

The clusters are saved in the field object$seurat_clusters where object represents your Seurat object.
You can also get them running the command Idents(object).
NOTE: if you run the command FindClusters on different SNN graphs, or different resolutions, you need to specify the clusters by object$graph_name_resolution.
Every time you run FindClusters, the idents of the cells are updated to the latest results, same for seurat_clusters.
For instance, let's say that you have run the command FindNeighbors on the RNA assay (default value).
Then you run FindClusters with graph = "RNA_snn" and resolution=1.
In that case, the clusters are in object$RNA_snn_1.
Once you have the clusters, you can save them as a CSV:
clusters <-  <here you use one of the described methods>
# data.frame
clust.df <- data.frame(cellID = colnames(object), clusters = clusters)
write.csv(clust.df, <path.of.the.csv.file>)

Thank you very much !


text:
How to convert amino acid sequence to SMILES format
I have an amino acids string "MKGK" and I would like to convert it into SMILES format.
I tried two approaches.
First with PepSMI which gives me:
N[C@@]([H])(CCSC)C(=O)N[C@@]([H])(CCCCN)C(=O)NCC(=O)N[C@@]([H])(CCCCN)C(=O)O

Second, using ChemAxon Marvin's molconvert command
/Applications/MarvinSuite/bin/molconvert --peptide MKGK smiles

I get:
[H]N[C@@H](CCSC)C(=O)N[C@@H](CCCCN)C(=O)NCC(=O)N[C@@H](CCCCN)C(O)=O

Both of them give different result. Are they chemically the same?
If not, which result should I trust?

SMILES are meant for machines, but reading up on the format is very handy and one totally can read and write them —I do it often and could not recommend it enough.
You can have SMILES written differently and they give the same result —hence why you should never compare SMILES, but Inchi or better still a parsed molecule —e.g. in RDKit.
Looking at them, I can see that

they are the same except for one explicit proton on the N-terminus.
Neither are canonical SMILES —which is fine.
Both are written N to C like biologist write peptides —chemists write peptides C to N due to solid phase synthesis.
they are L-amino acids: the pattern N[C@@H](*)C(=O)O is the give away. Technically, the implicit hydrogen in the Cα [C@@H] has @@, which means counter-clockwise  —Voet & Voet taught everyone to read C(=O)RN counter-clockwise with the N-C (the first bond) sticking out.

The first has rather oddly declared chiral protons, the second has an unnecessary N-terminus proton, which is not adding a charge: [H]N* is the same as [NH2]* as it is assumed to be [H][NH]*. So the correct version at pH 7 should be:
[NH3+][C@@H](CCSC)C(=O)N[C@@H](CCCCN)C(=O)NCC(=O)N[C@@H](CCCCN)C([O-])=O


text:
Build RepeatMasker default libraries beforehand (bioconda installation)
I'm using RepeatMasker on human sequences. When running RepeatMasker this way:
RepeatMasker -e rmblast -species human -dir . myFile.fa

RepeatMasker builds some libraries:
Building general libraries in: /path/to/share/RepeatMasker/Libraries/CONS-Dfam_3.3/general
Building species libraries in: /path/to/share/RepeatMasker/Libraries/CONS-Dfam_3.3/human

(Paths edited)
On subsequent runs, this does not happen.
I want to check if these libraries exist beforehand (which I can easily do with ls /path/to/share/RepeatMasker/Libraries/CONS-Dfam_3.3/), and if they don't, build them.
I could do a run with a decoy fasta to trigger that step, but it's a bit hacky. Is there a better way to force Repeat Masker to build those libraries?
(I mention that I'm using bioconda in the title because I'm working on a remote server with restrictions, so that should be taken into account, no because there is already a solution for other installation methods).

Something like this:
db_dir=$CONDA_PREFIX/share/RepeatMasker/Libraries/CONS-Dfam_3.3/

if [ -d "${db_dir}" ] && [ -n "$(ls -A ${db_dir})" ]; then
    echo "CONS-Dfam folder is not empty"
    # ... do something
else
    echo > sequence.fasta
    RepeatMasker -e rmblast -species human -dir . sequence.fasta
fi

If the folder exists and is not empty, otherwise build it using an empty file.

That's what I mentioned in my question, using a decoy FASTA. It seems there are not alternatives


text:
Blastn locally nt database search length < 16
when I run a short search in blastn locally with less than 16 nucleotides with nt.00 database I don't get any found match, but if I increase the searched string 2 nucleotides it finds matches, so the 14 string was in the BD. Also, when I run the same 14  nucleotides search in a BD created with homo-sapiens Y chromosome it finds matches. This BD was created with the command:
makeblastdb -in homoy.fasta -out homoydb -dbtype 'nucl' -parse_seqids -hash_index

As you see it includes -parse_seqids -hash_index parameters.
The command for search is (replacing homoydb with nt.00 in the first case):
blastn -task blastn-short -query myqueryin.fa -db homoydb -out myout.txt

So why blast nt preformatted databases don't find sequences shorter than 16 if they are in the BD and with my own DB find them OK?


text:
Generating PDB file for custom DNA sequence
I asked this question in biology.stack exchange and was referred to this community. I am going to request migration for this question.
https://biology.stackexchange.com/q/105024/67415
I want a PDB file of DNA sequence I am interested in (10 base pairs) for docking. I will appreciate if any of you can please clarify some of my doubts.
Do you think, I can generate PDB file by using this link http://www.scfbio-iitd.res.in/software/drugdesign/bdna.jsp
How can I get the pdb file of the DNA sequence I want from protein data bank, I did search it but did not find the one I want. Is it a way to extract the PDB file?
What are other methods I can use to generate PDB file of B-DNA sequence?
Thank you for your time.
Please delete if it is against the group for posting again and redirect me to where I can request migration of this question.

I normally use PyMOL. In PyMOL you can generate a model of a sequence of DNA with the command fnab — a sister command to the better known fab, which makes a protein sequence of a given secondary structure. I've not encountered that site, but it says it uses experimental data, which is a better approach.
Side note 1
OP is asking for DNA, but if anyone stumbles across this, do note that RNA is a totally different ballgame and requires more advanced solutions (example).
Side note 2
I should say docking a known sequence is fine, but not docking to guess the sequence is not.
For docking, if you have a zinc finger, helix-loop-helix, leucine zipper, MADS box or similar domain which binds DNA in a consistent manner then it is straightforward. However, predicting the bound sequence is far from trivial or accurate.
For example, in MEF2C the N-terminal amine there are three reader-head dictating specificity and one is a floppy terminus, which most forcefield-based methods would struggle to get right.
Testing
So let's say we want to test PyMOL's fnab in the GUI. We could pick a protein, MEF2A example from above:
fetch 1EGW, type=pdb1
print cmd.get_fastastr('chain E') # to get the sequence
fnab AAAGCTATTATTAGCTT
align obj01 and chain A, 1EGW and chain E
show sticks, chain E+F
hide cartoon
remove solvent

The RMSD is 1.6 Å which is not great, but pretty decent considering that it is with a generated sequence that does not take into account backbone distortion from the bound protein.

Thank you for your reply. I want to say that I know my sequence and the docking position. However, the big doubt in my mind is that if the pdb file generated through say pymol/avogadro/chimera is similar to the PDB file from the protein data bank? I want to learn how different it is if any.

Added a quick (nearly) coding-free test.

With probably the same results as PyMOL method there is x3dna webserver, creating an idealized DNA double-helix of customizable sequence: http://web.x3dna.org/fibermodel/regseq

PDB files from the PDB or generated are essentially the same in term of organization, though you would have a more complicated header and different B-factors for each atoms when the file comes from the data bank. Be careful not to mistake mmCIF format with PDB one. As for a difference in the structures Matteo answer covers it.

There is a tool for making DNA segments of any curvature. It exports as a PDB, and the people behind the tool are quite committed and helpful.
https://www.lifeexplorer.info/
Here is a curved DNA segment I generated using their tool; they can make much more dramatic bends.

I regularly use the site you cited (scfbio...) to generate idealized straight DNA segments for illustration purposes. It is not sequence-dependent. For an actual solved structure you can extract specific chains from existing structures containing DNA (e.g. RNA polymerase II structures by the Cramer group). These structures will have their own sequence, not your specific one, and will be distorted by protein-DNA interactions. Other than that, there likely won't be a PDB structure of your specific sequence. And even if there was, the conditions likely would differ from the ones your structure was acquired at. I may be entirely off, but it sounds like you should be fitting individual bases, sugars and phosphates, not continuous poly-d-nucleotides.
Shown below is the pymol visualization of a 100bp dsDNA input for B-DNA from the resource you cited It does not take into account the DNA flexibility (persistence length etc)


text:
uploading reads to ENA using aspera
Once again, I am trying to upload reads to ENA using aspera. Following their instructions I used several times in past (successfully) to upload my reads
ascp -QT -l300M -L- *.fastq.gz Webin-XXXXX@webin.ebi.ac.uk:.

it starts as usual, uploads by the ~300M/s speed. But when it gets to ~30 - 40% complete transfer, it gets aborted from the server site:
RR PEER ABORT: Received sesserr notification from peer, reason FASP_ABORT_LOCAL
ERR PEER ABORTED: Shutting down source/transmitter, reason FASP_ABORT_LOCAL
ERR Cannot write to management socket localhost:35666, Connection reset by peer (server)
LOG Sender bl t/o/r/dl/c/h=2672819/2672819/0/0/2668562/2668562 ctl bm/bs=0/0 rex b/r/n=0/0/0 r_ott q/s/n=0/0/1001 rate t/m/c=300000000/0/300000000 prog t/f/e=3762237952/376223795
2/104125014 xmtD=0
LOG Sender DS Qs ds/re/a/d/ao/eo/n/ap/rq/rd/ru/cl/di/do=1/0/0/0/0/0/0/0/0/0/1/0/0/0 Rs o/i=1/1 mgmt backlog i/s/n =
LOG Sender bl t/o/r/dl/c/h=2672946/2672946/0/0/2668562/2668562 ctl bm/bs=0/0 rex b/r/n=0/0/0 r_ott q/s/n=0/0/1001 rate t/m/c=300000000/0/300000000 prog t/f/e=3762237952/376223795
2/104125014 xmtD=0
LOG Sender DS Qs ds/re/a/d/ao/eo/n/ap/rq/rd/ru/cl/di/do=1/0/0/0/0/0/0/0/0/0/1/0/0/0 Rs o/i=1/1 mgmt backlog i/s/n =
ERR FASP sender had error, errcode 19 errstr [Cannot write to management socket localhost:35666, Connection reset by peer (server)]
ERR Peer FASP receiver had error, errcode 74 errstr [Cannot write to management socket localhost:35666, Connection reset by peer]
ERR Cannot write to management socket localhost:35666, Connection reset by peer (server), Cannot write to management socket localhost:35666, Connection reset by peer

Session Stop  (Error: Cannot write to management socket localhost:35666, Connection reset by peer (server), Cannot write to management socket localhost:35666, Connection reset by
 peer)

Does anyone have an idea what went wrong?
I also asked ENA helpdesk, if I will get an answer from them I will forward it here.

have you tried with a lower value for -l? Maybe -l200M

I tried -l100M and it failed around 15%, so I tried also 1500M and our bandwidth actually managed to transfer nearly 900M/s and although it dropped some of the same errors, it also said that the file transfer was successful. Might be a workaround, it's certainly not a scalable solution


text:
fastqCleaner failing to launch in RStudio
When I execute:
> FastqCleaner:::launch_fqc()()
I get the following output
Listening on http://127.0.0.1:4863
Warning: Navigation containers expect a collection of `bslib::nav()`/`shiny::tabPanel()`s and/or `bslib::nav_menu()`/`shiny::navbarMenu()`s. Consider using `header` or `footer` if you wish to place content above (or below) every panel's contents.
Warning: Error in : Navigation containers expect a collection of `bslib::nav()`/`shiny::tabPanel()`s and/or `bslib::nav_menu()`/`shiny::navbarMenu()`s. Consider using `header` or `footer` if you wish to place content above (or below) every panel's contents.
  85: stop
  84: FUN
  83: lapply
  82: buildTabset
  81: tabsetPanel_
  80: func
  78: tabsetPanel
Warning: Navigation containers expect a collection of `bslib::nav()`/`shiny::tabPanel()`s and/or `bslib::nav_menu()`/`shiny::navbarMenu()`s. Consider using `header` or `footer` if you wish to place content above (or below) every panel's contents.
Warning: Error in : Navigation containers expect a collection of `bslib::nav()`/`shiny::tabPanel()`s and/or `bslib::nav_menu()`/`shiny::navbarMenu()`s. Consider using `header` or `footer` if you wish to place content above (or below) every panel's contents.
  85: stop
  84: FUN
  83: lapply
  82: buildTabset
  81: tabsetPanel_
  80: func
  78: tabsetPanel

Have I done something wrong or is the function broken in some way?

The issue is with the navbar page, downgrading shiny to 1.6 or in my case 1.4 fixed it for me, the detailed link to find to the solution is here: https://github.com/rstudio/shiny/issues/3556
Cheers.


text:
Downloading exon-wise gene sequence using commandline
I want to download exon-wise gene sequence of around 200 genes from NCBI. Is there anyway to do this from terminal? I have the transcript id for each sequence and I am using efetch command to download their CDS.
eg: efetch -db nuccore -format fasta_cds_na -id NM_000125.4
However, I need their exon-wise sequence as well (i.e. I need the sequence of each exon separately). It will be better if the exon starts with start codon and ends with stop codon (the untranslated regions of those exons completely removed) but I can modify the sequences accordingly once I have all the sequences (UTR included).
Is there any way in which we can give the NCBI transcript id as an input and get exon-wise sequence as an output?


text:
The correct design for time series experiment
I have a treatment and control in two time points like this
> design
                                  condition time
    1       CTRL_24_hrs_replicate1   control   24
    2       CTRL_24_hrs_replicate2   control   24
    3       CTRL_24_hrs_replicate3   control   24
    4  treatment_24_hrs_replicate1        t   24
    5  treatment_24_hrs_replicate2        t   24
    6  treatment_24_hrs_replicate3        t   24
    7       CTRL_48_hrs_replicate1   control   48
    8       CTRL_48_hrs_replicate2   control   48
    9       CTRL_48_hrs_replicate3   control   48
    10 treatment_48_hrs_replicate1        t   48
    11 treatment_48_hrs_replicate2        t   48
    12 treatment_48_hrs_replicate3        t   48
    >

I want to test between treatment and control considering time point 24 hours to 48 hours
I have done like this
dds <- DESeqDataSetFromMatrix(countData=a,colData=design, design=~time + condition + time:condition)

But at this part I get error, although I am trying different things

> ddsTC <- DESeq(dds, test="LRT", reduced = ~ time + time:condition)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
Error in nbinomLRT(object, full = full, reduced = reduced, quiet = quiet,  :
  less than one degree of freedom, perhaps full and reduced models are not in the correct order

or
> ddsTC <- DESeq(dds, test="LRT", reduced = ~time:condition)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
Error in nbinomLRT(object, full = full, reduced = reduced, quiet = quiet,  :
  less than one degree of freedom, perhaps full and reduced models are not in the correct order
>

I tried these with no error although I am not certain if this makes sense at all while the results of both is the exactly the same
> ddsTC <- DESeq(dds, test="LRT", reduced = ~ time + condition)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
>

> ddsTC <- DESeq(dds, test="LRT", reduced = ~ time)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

Al I want is getting the difference of treatment versus control but considering time goes from 24 hours to 48 hours
Thanks for any help


text:
getting PDB file of interested sequence
I want a PDB file of DNA sequence I am interested in (10 base pairs) for docking. I will appreciate if any of you can please clarify some of my doubts.
Do you think, I can generate PDB file by using this link http://www.scfbio-iitd.res.in/software/drugdesign/bdna.jsp
I did generate the PDB file using the above site. However, when I am using that PDB file for docking it do not give accurate results. It just fits to the sequence and it lacks the information the PDB file downloaded from protein data bank has.
How can I get the pdb file of the DNA sequence I want from protein data bank, I did search it but did not find the one I want. Is it a way to extract the PDB file?
What are other methods I can use to generate PDB file of B-DNA sequence?
Thank you for your time.

Does this answer your question? Generating PDB file for custom DNA sequence

Partially not fully. I already did what was suggested by David.


text:
Identifying Variant Conservation Across Species
I've got gene sequences of mouse genes. Each sequence contains variants from the reference genome. I need to find out whether these SNPs/variations are in functionally significant regions (i.e. exons).
I also need to see if these SNPs/mutations are conserved across species.
I'm not very proficient in this, but would BLAST or ENSEMBL be useful here? Is there an easy way to search this? I have the sequences as .txt files. I'm proficient in Python if that's useful as well.

search for variant annotation tools. I can think of VEP in Ensembl, SnpEff, .. You would need your variants in VCF format (possibly other formats) and the tools annotate each variant with various qualitative/quantitative metrics you require

VEP online: https://www.ensembl.org/Tools/VEP


text:
How can I identify acceptor and donor atoms using BioPython?
I know how to traverse through structures, models, residues, and atoms in a protein chain using BioPython (which is very easy).
How can I identify donor and acceptor atoms in a chain?

donor/acceptor belongs to the main chain (MC) or side chain (SC) of amino acids ??

@pippo1980, side chain. My target is to measure the distance between donor and acceptor atoms.

Ligand or hydrogen atoms of other side chains ? Just trying to put it into a perspective

I guess once you know coordinates of each putative donor and acceptor, you calculate the distances between all the possible pairs and take the ones under a threshold. My PhD examiner was blabbing about right geometry too … but we were talking about X-ray crystallography no way to know where  the hydrogens are and a kind of middle resolution so … but remember geometry is important too

Here a informative link https://www.biostars.org/p/182325/ this links to another post : most software simply looks for donor + acceptor pairs separated by an appropriate distance and bond angle.

@pippo1980, How does DSSP/STRIDE help in finding donor/acceptors? They are just supposed to detect secondary structures, no?

Auch you got me !!! Used STRIDE to assign secondary on my models !!! Don’t know how pdf validate nowadays secondary structure signatures : here about stride http://webclu.bio.wzw.tum.de/stride/stride.pdf

DSSP. Dictionary of Secondary Structure of Proteins (DSSP) assigns eight state secondary structure using hydrogen bonds alone. DSSP defines a hydrogen bond where the bond energy is below -0.5 kcal/mol based on a Coulomb approximation of the hydrogen bond energy.

My hydrogen bonding problem was about Ligand —> side chain dono/acceptor. Now you can point to new way to asses that for ligand interaction too

A hydrogen bond donor heavy atom is one that is protonated, while the acceptor has a lone pair. If you talking about amino acids, then it's simple as they are fixed.

atom name
description
role
note

N
backbone nitrogen
donor
SS

CA
backbone Cα
apolar

C
backbone carboxyl C
sterically hindered and electrophilic
nucleophilic attack

O
backbone carboxyl O
acceptor
SS

CB etc.
sidechain carbons
apolar

OXT
C-terminal oxygen
acceptor
may be protonated especially if a membrane embedded C-terminus

Serine's OG / Threonine's OG1
sidechain hydroxyl
donor
SS: may form ST turns. Fairly uncommon: may be a donor

Asp OD1/OD2/Asn OD1
sidechain carboxyl/amide O
acceptor
SS: may form Asx turns

Glu OG1/OG2/Gln OG1
sidechain carboxyl/amide O
acceptor

Lysine NZ
sidechain amine
donor

Asn ND2 / Gln NG2
terminal admide N
donor
the direction of an Asn/Gln may be accidentally filled if manually assigned as N/O have the same density in X-ray crystals

Water can be both a donor and an acceptor.
Two charged residues form a salt bridge more so than a hydrogen bond, because the interaction is driven by the Coulomb electrostatic forces of the two heavy atoms.

A hydrogen bond donor heavy atom is one that is protonated, while the acceptor has a lone pair. -- Does this apply to BioPython's Atom object? I mean, can I test or search it using BioPython?

Yeah nah. A heavy atom is an atom that diffracts under X-ray crystallography —not hydrogens. The Atom instance is a graph network node representation of any atom. It may be a hydrogen —in silico protonated, neutron diffraction or 1H NMR. Technically to the most pedantically punctilious level, it could be virtual (e.g. Drude particle)...

So, can I use your supplied information and apply them on a PDB file via BioPython and extract donor and acceptor atom info from that PDB file?

Totally — assuming a pair is within about 3Å

So, if I can extract donor and acceptor atom-related info, can I predict Hydrogen bonds in this way?

Completely... as said somewhere else though it is way easier working in dihedral space and assigning the SS by Ramachandran plot bins, e.g. see ABEGO bins in SI Fig 1 of https://www.nature.com/articles/s41594-018-0141-6#MOESM1 —however, looking at the hydrogen bonds is more precise as you can better assign the correct type of turns.

@MatteoFerla so treshold for distance would be  about 3Å, does this change with resolution of the stucture ? what about the bond geometry (https://academic.oup.com/peds/article/15/5/359/1508259)?

@user366312 I guess you should add hydrogens to your pdb first

@pippo1980, I guess you should add hydrogens to your pdb first --- that option has been forbidden by my P.I. According to him, this adds too much noise to the generated data.

2.8Å is an ideal, >3.2Å is weak. Mathematically the Lenard-Jones term is used in thermodynamic simulations —in the wiki page you can see the potential plateaus quickly.  Sure, the true distance will be marginally affected by the b-factors of the atoms (both anisotropic/isotropic) but in solution the distance is likely a distribution (i.e. the atoms wobbling around).


text:
Linking UK biobank patient IDs with their corresponding genetic data
I am analyzing ukbb data and I have patient IDs for both cases and controls. Now I want the genetic data of all patients based on their IDs.
I am trying to test with chromosome 5 as I have both bgen and sample file. I am using plink 2 for the conversion of bgen and sample file
plink2 --bgen ukb_chr5_v3.bgen ref-first --sample ukb22828_c5_b0_v3_s487239.sample --make-bed --out testfile
I am getting three files as output like .pgen .psam .pvar.
This is how .psam looks like,

FID
IID
sex

482828
482282
1

and .pvar file has five columns, chromosome, position, ID, ref allele, alt allele.
I am still confused which IDs have the corresponding alleles/genotypes.
Please guide me in further analysis. My goal is find genetic data of specfic IDs which I am looking for.
Thanks

The genotypes are in the .pgen file, but encoded as a binary, so you can't query it directly. What are you going to do with the genotypes?

My aim is to calculate polygenic risk scores using pre-calculated weights for each individual and for that I possibly need effect allele count (0,1, or 2)


text:
In Pygraphviz, How can I assign edge colors according to edge weights?
I am now drawing a weighted directed graph using Pygraphviz.
The adjacency matrix A of the weighted directed graph is
A = pd.DataFrame({'X': [0.0,0.1,-0.8], 'Y': [-0.7,0.0,-0.1], 'Z':
    [0.4,0,-0.1]} ,index=["X", "Y", "Z"])

where A_{i,j} indicates the edge weight from the node i to the node j.
I was able to draw the graph without edge color.
import networkx as nx
from networkx.drawing.nx_agraph import to_agraph
from IPython.display import SVG, display

G_nx = nx.from_pandas_adjacency(A,create_using=nx.DiGraph())
G_nx.graph['edge'] = {'arrowsize': '1.0', 'splines': 'curved'}
G_nx.graph['graph'] = {'scale': '3'}

Agraph_eg = to_agraph(G_nx)
Agraph_eg .node_attr["height"] = 0.3
Agraph_eg .node_attr["width"] = 0.3
Agraph_eg .node_attr["shape"] = "circle"
Agraph_eg .node_attr["fixedsize"] = "true"
Agraph_eg .node_attr["fontsize"] = 8
Agraph_eg.layout(prog="neato")
Agraph_eg.draw('graph_eg1.png')

But I have no idea how to assign edge colors and widths according to the following two rules.
#1 If weights are more than 0, edge color should be red. If weight are less than 0, edge color should be blue.
#2 Edge widths should be proportional to the absolute values of weights.
Do you have any idea how to do this??

This might help you
for edge in G_nx.edges(data=True):
    color = "black"
    weight = edge[2]["weight"]
    if weight > 0:
        color = "red"
    elif weight < 0:
        color = "blue"
    edge[2]["color"] = color
    edge[2]["penwidth"] = abs(weight)

This gives you the following graph:

Depending on your network size this solution might take some time. "penwidth" is the graphviz information for the edge width.

Thank you veru much !! It works !!

@Apppii092 if my answer solves your question could you please mark it as accepted?


text:
What statistics can I use on qPCR data composed of several runs if all normality tests fail on the dataset?
I have >100 qPCR experiments that I'd like to analyze together, each containing the same set of genes (10 genes of interest and 2 reference genes).
I have four different samples (untreated & 3 treatments) from the same cell culture. Our setup allows two of them to be run at once, i.e. the samples from the same cell culture are split to two qPCR runs.

I would like to do statistics with the dCq values after normalizing to the reference genes, but the dataset doesn't pass any normality tests.
What methods are suitable for finding valid differences between treatments (or any other features like age)?

If the same amount of cDNA was used and reference genes were added to every run, is it still necessary to normalize between runs or between cell cultures?

A wilcoxon test is what you're looking for.
Yes, since there will inevitably be some small batch effects to compensate for.

Thank you! Could you give me tips about normalization because of batch effects?
Would you do a PCA as well to see which variables have the biggest effect on gene expression?

In qPCR the assumption is that there’s just a scaling difference between batches. So normally you load the same sample on multiple batches and scale according to that.


text:
Drug-DNA complex simulation using AMBER
I want to carry out the simulation of drug-dna complex by placing the drug molecule at particular site in my DNA sequence (intercalating site). I have PDB files of both drug and DNA sequence. I will truly appreciate if anyone could guide me how to place my drug at the desired location in the DNA sequence using Amber.

If it is only one... Can you not place it by dragging it in PyMOL?
Also do check out this discussion thread about the issue of aromatic dipole moments and forcefields: https://www.researchgate.net/post/Pi-pi-and-cation-pi-interactions-in-force-fields

Thank you so much for your reply. I can try that. However, I want to know how can I know  different base pairs in pymol. When I load the pdb file of my DNA sequence the base pairs cartoon all look same. Do you have any recommendations for teh software I can use to load my DNA and ligand pdb file and insert it into specific sequence and then save file as pdb  may be.

PyMOL shows everything as cartoon by default, but you can show them as licorice. Say show sticks, polymer.nucleic. There are a few other representations possible —cf. right hand side panel under (S) for show.

I should mention, to drag stuff, click on 3-button viewing on the right hand side panel making "editing", shift+wheel button ought to be Move object. Make sure to save often as there's no undo!

Thank you for your help


text:
Constructing a phylogenetic tree for classical MHC genes from different organisms
I am interested in constructing a phylogenetic tree containing the classical HLA α1 and α2 domains and potential classical MHC α1 and α2 domains of a different organism. I am using the amino acid sequences for the tree.
However, I am unsure of what bioinformatic tools and software to use and what parameters to fill in.
Lastly any advice on how to read a phylogenetic tree and how to relate the sequence similarity between two sequences by looking at the tree would be very helpful.
My aim is to investigate the classical MHC genes of a non human organism.
Any help and advice would be greatly appreciated!

This is a very unusual situation as a main effort is finding the outgroup where these two domains aren't duplicated.
General situation
In animals a gene tree nearly always follows a species tree, with the exception of gene duplication events, which frequently stem from genome duplication events. Making a tree with lots of genes gets hard to follow and actually becomes pointlessly problematic when you get a branch that does not quite follow the species phylogram —say, you have rabbit, mouse and human genes and the inferred tree, due to poor signal, comes out with mouse+human together as opposed to rabbit+mouse. As a result a nice tree needs to present the species before the genome duplication, often sea lamprey —I wrote a blog post recently were I discuss nice species to pick for a multiple sequence alignment.
In this simple case, one can get the protein sequence of interest and blastp it against human, sea lamprey and fruit fly, see how divergent are the paralogues in human and if the split is more recent that the latter two, if so, the latter will be an outgroup and after getting a few representative species in between, one can make a MSA with Muscle, which contains also a tool to make a neighbour joining tree.
Domain duplication
In your case you have something more complicated in that you have a series of paralogues, each gene encoding a protein (e.g. HLA-A) containing two globular domains, which form the MHC-1 fold, and you want make a tree of this... which means you want find out when the two domains duplicated (i.e. outgroup)? If the domains work together (correct?) and have the same fold (as you alluded to, correct?), this would mean that before the duplication event of the domains, the ancestral protein may have worked as a homodimer.
I am not familiar with the literature, so do not know if this has already been looked at —but I suspect it has.
If not, clicking on the sequence of one of the two domains in the Uniprot entry of one of the paralogues and doing a few iterations of psi-blast against human and a few other distant species (such as zebra fish, sea urchins, which are basal to vertebrates, and fruit fly, which I do not believe have a MHC etc.) most of these will have two hits per protein... but one might not. This is the outgroup!


text:
Using python to call Emboss wordcount
Hi guys I would like to run wordcount from Emboss in a bunch of fasta files. I tried this:
#!/usr/bin/env python
import sys
import subprocess

filenames = ['Downloads/Data/GCF_000016725.1_chr.fna']
for filename in filenames:
    name = filename.split('/')[2]
    print(filename)
    subprocess.call (['wordcount', 'Input sequence(s): ' + filename, 'Word size [4]: ' + sys.argv[1]])

But I got error messages:
Count and extract unique words in molecular sequence(s)
Error: Failed to open filename 'Input sequence(s): Downloads/Data/GCF_000016725.1_chr.fna'
Error: Unable to read sequence 'Input sequence(s): Downloads/Data/GCF_000016725.1_chr.fna'
Died: wordcount terminated: Bad value for '-sequence' and no prompt

But if a use shell with Downloads/Data/GCF_000016725.1_chr.fna  it works fine!
Any tip would be helpful!
I am using python, but I would have a method to compare my results!
Thank you by your time!
Paulo

You have to use command line options, not interactive prompts, like so:
subprocess.call(['wordcount', '-sequence', filename, '-wordsize', sys.argv[1]], '-outfile', filename + '.wordcount.txt')

Or, using the f-strings or formatted string literals:
subprocess.call(f'wordcount -sequence {filename} -wordsize {sys.argv[1]} -outfile {filename}.wordcount.txt')

Note that I had to add -outfile option to your system call, otherwise wordcount prompts for the output file. See wordcount manual.

I tested your first solution and worked very nicelly. Thank you by your time and patience. Paulo


text:
How to show electrostatic interactions in Pymol on the wanted residue
I have the following Pymol visual,
The green molecule is receptor and red/yellow is the ligand.
In the ligand I highlighted Lysine (K) residue as yellow.

What I want to do is to highlight and draw the electro static contact of the ligand especially with Lysine residue. Like this in Pymol:

How can I achieve that?

Depends whether you work with named selections which is preferable and helpful if you want to use this Pymol session for future images. In this case, use the distance command for your two selections, one being the N of the K in your ligand and the other the O of the E of your protein.
Full syntax and examples here: https://pymolwiki.org/index.php/Distance
If you need a quick solution: in the menu Wizard, click on Measurement. Then click on the two atoms where you want the dotted line to appear in between.

Thanks but I'd like to see the electrostatic contact - proton and electron - like shown in my OP.

This answer is still valid, you simply click the hydrogen not the nitrogen for your distance. However, I would suggest that you do color atomic, not element C as you currently have all elements coloured in and it's well confusing...


text:
Can you style network visualizations of directed graphs using NETWULF?
I want to draw directed graphs using networkx, followed by manual tuning of the visualization with the use of NETWOLF.
Here is the example of a network.
G =nx.DiGraph()
G.add_edges_from([(1, 2), (1, 3),(1, 4),(2, 3)])
network, config = nw.visualize(G,config={'zoom':3})

However, the graph shown in NETWOLF was undirected graph.

Do you know if you can style network visualizations of directed graphs using NETWULF?

I'm the author of gravis, an open-source package for interactive graph visualization in Python. It allows to draw directed graphs from NetworkX and the resulting visualization can be modified by user interaction (e.g. node positions, edge curvature and various other visual properties).
Here's a minimal example based on your code:
import gravis as gv
import networkx as nx

g = nx.DiGraph()
g.add_edges_from([(1, 2), (1, 3), (1, 4), (2, 3)])
gv.d3(g)

Output: NetworkX DiGraph visualized with gravis inside a Jupyter notebook:


text:
I have two fasta files READ1 and READ 2. How to identify which one is the reverse and the forward?

As you can see in the picture I have two files from the paired end reads. I would assume READ1 would be forward reads and READ2 would be reverse reads, but my lab colleague told me that that might not be the case. How do I confirm?

Hi @fullmooninu - any idea what type of data these reads are? It will help us give you an accurate answer.

Unless this is RNA-seq prepared with a directional library, or something else that you're sure is directional, it is irrelevant.
DNA is complementary, so a DNA molecule can be turned into a sequence-able fragment with either strand given any pair of Illumina adapters and indices.
As far as Illumina lingo goes, R1 is the first read that is read off the sequencer, or what people call the "forward read". R2 is the second pass of the same clusters, and some people call it the "reverse read". I recommend just calling them R1 and R2 reads to avoid confusion. Sequencing by synthesis produces the cleanest reads in the first pass, so R1 reads usually have better Q values. You can check this easily by comparing the quality scores with fastqc.
You can learn more about Illumina sequencing chemistry here, or check out some other answers on the bioinformatics SE or biostars.

I think @conchoecia answer is correct but maybe not clear enough for someone making himself that question for the first time (like was my case recently).
I think the simple way to explain this is that in F_NL_READ1.fq you have random reads, in mixed directions, and in F_NL_READ2.fq you have the paired read in the opposed direction for every read in F_NL_READ1.fq, but not necessarily and almost certainly, not every read in every file will go in the same sense (forward or reverse).
At least for paired-end reads from Illumina.

Since we have no idea what kind of reads the OP has, what you wrote might very well be completely wrong.

If this is RNASeq, then the question might have an answer (or it might be nonsense), but you have to know what kind of library prep was used.
Or just align and inspect the results to figure it out.