text:
BWA: Detecting Variation between Reference Genome and Short-Read Sequences
I need to identify all loci in the short-read sequence at which the number of microsatellite repeats (i.e. number of copies of "AA," "GTC," etc.) differ from the reference genome, as well as the locations at which the number of repeats is the same in both the reference genome and the short-read sequence.
I used the Burrows-Wheeler Aligner (BWA mem) to map a high-coverage short-read sequence (obtained from the NCBI's short-read sequence archive) to a reference genome. The output is in .sam format. I have also used a separate program to identify the loci in the reference genome at which microsatellites occur.
I would like to identify all the loci in the short-read sequence at which the microsatellite length and loci differ from the reference genome. Does anyone know any tools or packages I could use to read a .sam/.bam file of a short-read sequence mapped to a reference genome and identify specific loci at which the short-read sequence differs from the reference genome? I am using RStudio and have access to my university's supercomputer clusters.
For info on microsatellites, see here: https://en.wikipedia.org/wiki/Microsatellite#:~:text=A%20microsatellite%20is%20a%20tract,locations%20within%20an%20organism's%20genome.
My end goal is to determine the loci at which microsatellite alleles tend to differ within a species, as well as the loci at which microsatellite alleles tend to be identical between all members of the same species.

Can you please provide more information about your specific problem? This will greatly help in working through a solution. There is presumably some end-goal that this STR detection is used for. No algorithm will detect all variants, most will produce some amount of false negative and false positive results.

I haven't done any whole-genome STR analysis from NGS data myself, but are aware of others that have used lobSTR for this. There's also a recent paper [here] that compares a few different STR analysis packages (i.e. RepeatSeq, LobSTR, HipSTR, GangSTR). Here's the concluding paragraph:

In conclusion, all these tools are built to genotype STRs but
have different strengths and weaknesses. Based on our analysis there is
no clear overall winner. RepeatSeq and HipSTR are the best when
considering genotyping error rate even with low coverage. On the other
hand, GangSTR has an advantage because it is the only tool among them
that can call alleles longer than the read length but shows higher error
rate unless looking at only the enclosed class of reads, which in turn
would lose the GangSTR's advantage of picking up long genotypes. In
addition, GangSTR is the newest tool and so comes with reference files for
different reference builds that are periodically updated according to
the tool's webpage. The correct choice of a tool and the
subsequent filtering depends on the aim of the analysis, and might be
influenced by available hardware resources and time limit for running
tools.

In other words, without more information about the specific problem you're trying to solve, it's not clear what the best tool to use is (this is a common issue with Bioinformatics problems).

ah I see - thank you so much for linking the paper! I will give it a read. My specific problem is that I need to identify all loci in the short-read sequence at which the number of microsatellite repeats (i.e. number of copies of "AA," "GTC," etc.) differ from the reference genome, as well as the locations at which the number of repeats is the same in both the reference genome and the short-read sequence

Honestly, if anyone knows of a software that works on R and detects indels using BWA-generated files, that would work for my purposes

Yep its a cool answer, I like it too. Your comment above is a separate question. My personal thoughts are you'd hit a RAM bottleneck pretty quickly via R and you should think about a more complex data pipeline.


text:
Is there a way to customize what observations are plotted in ComplexHeatmap?
So right now my observations in the dataset I've processed is stored as sample identifier columns and gene rows: an example of this would be for sample 1 and gene TRIM21 the observation is Missense Mutation;Frame Shift Deletion. I'm having issues separating the two alterations and instead of being able to show the gene in the sample having 2 different alterations, I end up having to generate a new colour classification (being that the gene in the specific sample exhibits specifically both a missense mutation and a frame shift deletion instead of being able to represent the observation with just colours and missense and frame shift deletion colours alone). Here is an a screenshot of what my matrix looks like:
As an example here is some test code generating the test plot:
library(ComplexHeatmap)

test_labels <- c("Missense_Mutation", "Splice_Site", "Frame_Shift_Del",
                 "Splice_Region","Missense_Mutation;Splice_Site", "Splice_Site;Frame_Shift_Del",
                 "Splice_Site;Missense_Mutation", "Frame_Shift_Del;Splice_Site")

plot <- Heatmap(test, row_names_side = "left",
        row_names_gp = gpar(cex = 0.7),
        column_names_gp = gpar(cex = 0.7),
        heatmap_legend_param = list(title = "Alterations",
                                    at = test_labels, labels = test_labels), show_heatmap_legend = TRUE)

plot

I've attached a screenshot of my actual plot and as you will be able to see, this ends up creating multiple different classifications:

Dear John, I am sure your post makes sense to someone, but makes little sense to me. I can see in your tags you are using Bioconductor ... that is a good start. Perhaps you can post your code and explain the data set. I see that the red is stuff you don' twant us to see

Hi Michael, thanks for responding! This is specific to a package in Bioconductor: https://jokergoo.github.io/ComplexHeatmap-reference/book/index.html. I would add my code, but I would also need to provide the exact dataset as well as the data processing I've already done, which is a decently long script. Essentially, the data set is represents samples (columns) and genes (rows), and observations correspond to what the alteration present in the specific gene in the specific patient sample. As you can see, sometimes observations are separated by ";".

This is a pretty specific question I realize, as I am using the Heatmap function in the vignette and it doesn't exactly cover how I'd be able to parse what I'd like to.

Thanks. Don't need to provide your data set ... a dumbie data set is fine - we understand patient confidentiality issues. Putting some code there would be useful.

I've updated my question with some sample code! You'll see in the legend, e.g. "Missense_Mutation;Splice_Site" is given its own category instead of being recognized separately as "Missense_Mutation" and "Splice_Site".

What color do you expect when the algorithm encounters a value that is neither Missense nor Frame-shift but Missense;Frame-shift? Are you looking for the cell to be split into 2 colors if that happens?

Because if that's the case, eyeballing will be extremely difficult with each cell looking a different height (or width) and the plot showing jarring patterns. Pick the most severe consequence and plot those alone, or use cell_fun to create a plot that IMO would be a mess and beats all reason for visualization.

What is fed into ComplexHeatmap is what it outputs! So if one of your genes has a Missense_Mutation as well as a Frame_Shift_Del and you would like to plot its value, of course you will have a category and color that is different than those of Missense_Mutation or Frame_Shift_Del, that is only natural.
What you can do is to create additional variables like your_gene_Missense_Mutation or your_gene_Frame_Shift_Del each with 0 or 1 to indicate if a particular mutation type is present or not for a particular gene (you can think of it like one-hot-encoding). With a lot of different levels, the resulting heatmap could be difficult to interpret though.
Another simpler alternative would be a using a "universal label" if your gene has multiple mutations. Basically if the cell contains a ;, its value will be changed to something like mixed_mutation_status and instead of multiple "mutation_type_1; mutation_type_2; ... like categories, you will just have one mixed_mutation_status on top of Missense_Mutation, Frame_Shift_Del, ...

The mixed_mutation_status is indeed how maftools does it - it calls it multi-hit. Although, in this case, the same mutation can be both missense and splice_site if it's a mutation at an exonic spllice site. OP should first pick among the many possible consequences of a mutation wherever possible, and then apply a logic where multiple mutations still exist in the same gene.

source: Stack Exchange

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

sort: votes


text:
Remove Redundant Sequences from FASTA file in Python
I'm attempting to remove redundant sequences from a fasta file (from NCBI). When I execute this code, it returns the number of spots, not the number of sequences. (Number of spots: 408,293, Number of sequences: 40,830)
Please let me know if you have run a code that can remove the duplicate ideas or have any ideas/suggestions.
Thank you SO much.
from Bio import SeqIO
import time

start = time.time()

seen = []
records = []

for record in SeqIO.parse("b4r2.fasta", "fasta"):
    if str(record.seq) not in seen:
        seen.append(str(record.seq))
        records.append(record)

#writing to a fasta file
SeqIO.write(records, "no_dupes_b4r2.fasta", "fasta")
end = time.time()

print(f"Run time is {(end- start)/60}")

Hi, your code was not formatting and it became a bit of a headache to format. Please invest some more effort into presenting your questions better.

Thank you so much, I'm just learning to use Stack Exchange! I appreciate it

Your code (as currently formatted) should generate the correct output, but it is not efficient. Make seen a set to improve these O(n) in calls to O(1) time complexity- this will speed it up a lot
seen = set()
records = []

for record in SeqIO.parse("b4r2.fasta", "fasta"):
    if record.seq not in seen:
        seen.add(record.seq)
        records.append(record)

Thank you for your input!

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


text:
What is the origin of the 'source Gene ID' references given in the 'gene_presence_absence.csv' output of Roary?
I am learning to use Roary for preparing a pan genome for some lactobacillus strains.  In the 'gene_presence_absence.csv' output of Roary (which I view in excel), a 'source Gene ID' is given for each of the strains included in the analysis.  The Roary web page states that this is:

Presence and absence of genes in each sample, with the corresponding source Gene ID.

Is this 'source Gene ID' generated by Roary, or does it relate back directly to the .gff files I used as input?  If I take one of these 'source Gene ID' references and Grep in the directory containing my .gff files I get no matches, so I am assuming it is generated by Roary...?
I had hoped that this reference would link back directly to the gene within the annotation file, but that seems not to be the case?  Should/can I use this reference for anything?

Did you end up getting an answer to this question? I have the same question.

@Sharb - I think i can now answer this question.  I have added an answer.  Let me know if you fee it needs any clarification?

Looking at this again, I have now worked out (and should probably have noticed previously) that the 'source Gene ID' field from the 'gene_presence_absence.csv' file originates from the ID field of the .gff file for the sequence in question.
You can therefore easily tie the results from the 'gene_presence_absence.csv' file back to the annotation from the sequence in question, and therefore to the location of the nucleotide sequence for the gene in question.
So for instance, if you take a 'source Gene ID' from the gene_presence_absence file - say GEICCBEJ_00950 - and use it in a grep command referencing the genome from which it came, grep "GEICCBEJ_00950" my_annotation_file.gff, then you should see the line from the annotation corresponding to this sequence.
The 'source Gene ID' is also referenced in the pan_genome_reference.fa output from ROARY, which gives "single representative nucleotide sequence from each of the clusters in the pan genome (core and accessory)".


text:
What does aov(glm model) actually do in R?
What is happening when you use aov() on a glm model in R?
Normal glm model sumamry(m1):
Call:
glm(formula = value ~ origin + variable, family = poisson(),
    data = comb_data)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-0.8227  -0.5547  -0.4385  -0.3038   2.5258

Coefficients: (1 not defined because of singularities)
                                                                  Estimate Std. Error z value Pr(>|z|)
(Intercept)                                                     -2.853e+00  5.774e-01  -4.941 7.78e-07 ***
originisolate                                                    4.700e-01  7.071e-01   0.665  0.50625
variableRB17                                                     2.877e-01  5.401e-01   0.533  0.59425
variableM56                                                      1.299e+00  4.606e-01   2.821  0.00479 **
variableRB68                                                    -6.269e-14  5.774e-01   0.000  1.00000
variableM39b4                                                   -1.010e-13  5.774e-01   0.000  1.00000
variableRB108                                                    1.542e-01  5.563e-01   0.277  0.78172
variableM32                                                      4.055e-01  5.270e-01   0.769  0.44171
variableRB33                                                    -6.931e-01  7.071e-01  -0.980  0.32696
variableM41                                                      6.061e-01  5.075e-01   1.194  0.23236
variableRB1                                                      6.931e-01  5.000e-01   1.386  0.16566
variableM39w2                                                    1.542e-01  5.563e-01   0.277  0.78172
variableRB20                                                    -1.792e+00  1.080e+00  -1.659  0.09714 .
variableRB99                                                     1.542e-01  5.563e-01   0.277  0.78172
variableRB110                                                    6.931e-01  5.000e-01   1.386  0.16566
variableRB5                                                      1.542e-01  5.563e-01   0.277  0.78172
variableRB29                                                    -1.823e-01  6.055e-01  -0.301  0.76334
variableRB56                                                    -7.411e-14  5.774e-01   0.000  1.00000
variableRB24                                                    -1.099e+00  8.165e-01  -1.346  0.17846
variableStreptomyces_alboflavus_strain_MDJK44_RB5                1.204e+00  6.583e-01   1.829  0.06740 .
variableStreptomyces_violaceoruber_RB1                           1.299e+00  6.513e-01   1.995  0.04607 *
variableStreptomyces_T44_RB17                                    1.204e+00  6.583e-01   1.829  0.06740 .
variableStreptomyces_violaceoruber_RB1_RB110                     1.299e+00  6.513e-01   1.995  0.04607 *
variableActinomadura_verrucosospora_RB99                         5.108e-01  7.303e-01   0.699  0.48425
variableStreptomyces_autolyticus_strain_CGMCC0516_M56            1.735e+00  6.262e-01   2.770  0.00561 **
variableStreptomyces_atratus_RB108                               9.808e-01  6.770e-01   1.449  0.14740
variableStreptomyces_sp_T44_RB17                                 1.204e+00  6.583e-01   1.829  0.06740 .
variableStreptomyces_sp_QMT.28_M41                               8.473e-01  6.901e-01   1.228  0.21950
variableAmycolatopsis_sp_AA4_M39                                 8.473e-01  6.901e-01   1.228  0.21950
variableLuteimicrobium_xylanilyticum_RB22                       -1.099e+00  1.155e+00  -0.951  0.34138
variableMycolicibacterium_parafortuitum_RB24                     1.239e-13  8.165e-01   0.000  1.00000
variableNocardia_nova_SH22a_RB20_RB56                            2.877e-01  7.638e-01   0.377  0.70642
variableActinomadura_atramentaria_RB29                           5.108e-01  7.303e-01   0.699  0.48425
variableMycolicibacterium_madagascariense_strain_JCM_13574_RB33         NA         NA      NA       NA
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 1024.95  on 1949  degrees of freedom
Residual deviance:  943.36  on 1917  degrees of freedom
AIC: 1507.4

Number of Fisher Scoring iterations: 6

A summary(aov(m1)):
              Df Sum Sq Mean Sq F value   Pr(>F)
origin         1   0.33  0.3285   3.050   0.0809 .
variable      31  10.41  0.3357   3.117 2.06e-08 ***
Residuals   1917 206.47  0.1077

I’m voting to close this question because 1) the OP only throws in a R function and tables without a reproducible example, 2) it seems to be mostly about statistics and not anything to indicate it is about bioinformatics

welcome to bioinformatics SE. If you indeed have a question related to bioinformatics, please provide more context and what you would like to do. See also GordonSmyth's answer below. You normally compare deviance to see how useful are the terms, not the sum of squares as you would for a linear model

I am attempting to gain an understanding of that is happening inside an R function, Cross Validate don't like R specific questions so I came here. I would argue this is a bioinformatics related questions and a reproducible example is not needed to answer the question I pose. It would just further clutter the question. However I guess some context could have helped.

It's not appropriate to run aov on a GLM model. aov assumes normally distributed data and it runs on a GLM as if it was a least squares linear model. It uses the working responses and working weights from the GLM fit and performs an ordinary anova as if the responses were normally distributed.
On the other hand, anova is a generic function with a method for GLMs and it performs a sequential analysis of deviance.
Purely statistical questions like this are best sent to CrossValidated, although they don't like questions that are posed entirely in terms of an R function. Questions purely about the functionality of functions in the core distribution of R are best sent to R-help, which is a mailing list rather than web forum (see https://www.r-project.org/mail.html and https://stat.ethz.ch/mailman/listinfo/r-help).

It's equivalent to doing anova(m1, test="F")
You probably want to do some sort of deviance test for a Poisson glm
anova(m1, test="Chisq")

But what does performing an anova on a glm model do, and why does this one throw out a p value for each parameter of the model whilst anova() does not?

Probably best to ask on CrossValidated. If you use the syntax in my answer, you'll see that p values are given with anova. Again if you want to know why anova doesn't give p values by default for a glm best to ask on CrossValidated - or search anyway...


text:
e-utils api: get chromosome position for each snp on a gene
So I have the following e-utils api-link:

https://eutils.ncbi.nlm.nih.gov/entrez/eutils/elink.fcgi?dbfrom=gene&db=snp&id=5726&retmode=json

it returns all the SNP-ids for the gene with the id 5726 as an example. Additionally to the SNP-ids I also need the SNP position and the genomic distance. Also I would need the chromosome but I could query that from the gene as well rather than for each SNP. Do you know how I can add all these informations to the resulting json/xml?
The first one in the resulting list is the SNP with the id 1490381980. Further information about it can be found on dbSNP: https://www.ncbi.nlm.nih.gov/snp/?term=1490381980%5BUID%5D
I really want to avoid creating a request for each single SNP and would prefer to query everything at once. Due to the unavoidable performance-issues. The Tas2R38 doesn't have that many SNPs but there are other genes with way more.

Hi @cowboy_patrick could you plase state/clarify what you mean by "to query everything at once. Due to unavoidable performance issues.".

It would be great if I could get all the SNPs and their positions with a single api-call for a specified gene. Instead of creating a call for each SNP on the variation services or something similar

If your dataset is limited to human SNPs, and you know the list of SNPs, I would advise taking a look at the dedicated NCBI Variation Services.  This would be one request for each SNP, but if you fire off 1-2 per second, I think you'll be ok.
For example, for rs1490381980, you could try something like:
curl -X GET 'https://api.ncbi.nlm.nih.gov/variation/v0/refsnp/1490381980' \
   | jq '.primary_snapshot_data.placements_with_allele[] | "\(.seq_id) \(.alleles[0].allele.spdi.position)"'

It would be better to query for the SNP and parse the JSON in python, but jq helps provide an idea of what is possible.
Alternatively, you could download the JSON or VCF data from the NCBI dbSNP FTP site, and process out the information you need.


text:
Calculating average coverage for .bam files (sequence data)
(Full discolosure that this is my first time working with sequence data, and with the bash scripting.)
I need to calculate the average coverage for any .bam file.
After some searching I wrote the following script:
# Script to calculate the average coverage of a genome sample

# compute the total length of the sample
tot=$(samtools view -H $1 | grep -P '^@SQ' | cut -f 3 -d ':' | awk '{sum+=$1} END {print sum}')
echo $tot

#compute the coverage at each point
samtools depth -ao cov_$1 $1 #your .bam file has depth calculated for each point

#calculate the average coverage
sum=$(awk '{sum+=$3} END {print sum}' cov_$1)
echo $sum
echo

avg=$(echo "$sum/$tot" | bc -l)
echo ''The average coverage of the sample $1 is $avg x.''

The logic behind the calculation of tot is as follows:
samtools view -H  gets the header of the bam file
grep -P gets all lines in this header with SQ
cut -f 3 -d ':' removes fields that are lists of more than 3 letters while -d uses delimiter ':' instead of tab and this leaves only the numbers i.e. lengths of the chromosomes
awk '{sum+=$1} END {tot=sum}' calculates the length/total bases in the chromosome
However, my output seems incorrect. Any advice on where I might be making an error (in code or thought process), or a better method would be appreciated.
One thing that worries me is that my tot variable is the total number of bases, so if the samtools depth function does not calculate coverage for each base then my sum and tot will not be in the same units and my calculated coverage will be wrong...

For a quick estimate you’re making it more complicated than necessary.

The theoretical average coverage is $\frac{n \cdot \hat l}{N}$ where $n$ is the number of reads, $\hat l$ is the average read length and $N$ is the genome size.

samtools idxstats gives you the chromosome lengths and number of mapped reads in one convenient list.

Putting this together, we get
samtools idxstats data.bam \
| awk -vreadlen=120 '
    {
        len += $2
        nreads += $3
    }
    END {
        print nreads * readlen / len
    }
'

(In this case I’ve used 120 as the read length.)
However, note that this over-estimates empirical sequencing depth, mainly due to the fact that a given read may not be fully mapped; thus the the average read length is an upper bound on the average mapped length. If you have an estimate for the average percentage of mapped bases per read, you can multiple that by the above coverage estimate for a more accurate number.

If you need a more accurate estimate, your approach is fundamentally sound (and the result is probably correct!). For low-coverage data you can speed up your analysis by dropping the -a argument to samtools depth: you don’t need to record zero-coverage bases, they don’t contribute to the total. You also don’t need to store the coverage into an intermediate file, thus reducing IO cost:
sum=$(samtools depth "$1" | awk '{sum+=$3} END {print sum}')

Some more remarks on your code:

The -P flag for grep is not supported everywhere and thus limits portability, so you want to avoid using it if possible. It turns out that you don’t need it here: grep '^@SQ' works just as well.
While the philosophy of UNIX is based on solving problems by combining individual tools that each solve one part, it is sometimes convenient to merge individual steps:

You can often combine cut and awk, in this case by setting the field separator:
awk -vFS=: '{sum+=$3} END {print sum}'

You can often combine grep and awk:
awk -vFS=: '/^@SQ/ {sum+=$3} END {print sum}'

Always double-quote your shell variables, otherwise your script will break if the values contain unexpected spaces or other special characters. This is especially important for filenames, because spaces and other special characters are valid characters in paths!
Your code is clean and readable. Well done.

With the current versions of samtools, you can also use:
samtools coverage sample.bam

The documents are here,

samtools coverage – produces a histogram or table of coverage per
chromosome

A further example from the site,
samtools coverage -r chr1:1M-12M input.bam

Output
#rname  startpos  endpos    numreads  covbases  coverage  meandepth  meanbaseq  meanmapq
chr1    1000000   12000000  528695    1069995   9.72723   3.50281    34.4       55.8

The histogram option is
samtools coverage -A -w 32 -r chr1:1M-12M input.bam

... and output available on the website link.


text:
Bash script for automatic overlap merging of respective paired-end reads from samples
I have around 300 paired-end files from 150 samples (1 forward and 1 reverse read for each sample). I want to merge (by sequence overlap) respective forward and reverse reads for each of the samples. The samples are named e.g. SRR2155174_1.fastq SRR2155174_2.fastq SRR2155319_1.fastq SRR2155319_2.fastq. Can anyone please help me write a bash script to do this?

Hi @CHANDRA7, we do understand "paired ends", we'd like to know reaad-length, coverage, genome - don't have to but it is polite. HOWEVER, whilst there is serious level coding skill on the site fior e.g. bash/awk and its not hard (even I code this one)... this is not a code writing service and you will need to at least make an attempt at a line of code .I suspect you want to concatenate e.g. `SRR2155174_1.fastq SRR2155174_2.fastq > SRR2155174_12.fastq'. Please at least try (hint a bash do loop) - someone can code the regex part.

Cross-posted: https://www.biostars.org/p/447550/
Michael is telling you basically the same as I did over at biostars.

Since you have many files it may be helpful to use tools like snakemake to run automatically, although as previous comments said the question has been answered elsewhere.

Why do you want to concatenate them? That's almost never the right thing to do.

Respected @DevonRyan , I am doing microbial profiling of metagenomic WGS sequence with metaPhlan software which can't generate contigs. But they can do profiling while we catenate the forward and reverse reads. Due to its short marker sequences it can do profiling from unpaired reads and probably average out. Please check this: https://groups.google.com/forum/#!topic/metaphlan-users/aclhwUcgWkk . And please let me know if I have wrong understanding.

Your original question mentioned concatenating things, which is not what you want (that would result in all of the reads from the _2 files being after all of those in _1). What you want is to merge overlapping paired end reads. There are a variety of tools available for that, from flash to bbmerge. The general structure then is:
for R1 in *_1.fastq ; do
    sampleName=${R1%%_1.fastq}
    some command with ${sampleName}_1.fastq and ${sampleName}_2.fastq producing ${sampleName}.fastq
done

You can figure out the some command ... portion yourself.


text:
how to have a blast xml file in a readable and understandable way
I got a blast xml file after a DNA sequence search. The XML file is not very readable. What is the way to make it more understandable? This blast contains about twenty matches. Here is the xml file

You can do this using the Bio Entrez parser in Biopython. Its quite easy to use, but its not clear whether your happy using Python.

The most human-readable blast output formats are 0-4, e.g. -outfmt 0 - which if I am not mistaken is the format used by web blast. However, these formats are a pain to automatically parse.
As Cowboy_Patrick pointed out, XML is more of a machine-readable format. The tabular output format ( -outfmt 6) is very commonly used, because it is easily parsed by unix command-line tools, and its structure also allows a somewhat easy human reading and understanding.

Many thanks, how to set the output format as -outfmt 6 in a line of python code as
result = NCBIWWW.qblast('blastn','nt',ADN.format('fasta') ) ?

@serpipanton I am unfamiliar with BioPython, but a cursory look at the documentation suggests -outfmt 6 is unavailable for NCBIWWW.qblast. However, there is a module NCBIXML, which will allow you to easily parse the xml output.

Is this what you are looking for: https://codebeautify.org/xmlviewer ? XML is a format designed to be readable firstly by machines, not humans.
With xslt you can convert XML to various more readable formats see this blog for something that fits your XML: http://plindenbaum.blogspot.com/2008/05/ncbi-blast-xslt-xhtml-svg.html?m=1


text:
Mapping UniProt id to Entrez id
How can we get mapping from UniProt id to Entrez ID? I have UniProt id, but want to convert those ids to Entrez or Ensembl id.

The FTP sites of UniProt and Ensembl have mapping tables.

If you use R and you have a model organism, there are bioconductor packages like org.Hs.eg.db that allow you to map

If you can reformulate the question @user224050 to include why this is relevant to you and your investigation. There are perfectly good answers here, but right now the question is likely to 'bomb out', most likely get closed because anyone who has not encounteed the problem would not understand the question and it must be understandable to a broad audience.

I think this is exactly the type of question we need - concise and googleable. I just took the liberty of making the title more informative (hope OP won't mind).

Baring the response being "accepted" for the answer below, I agree the question worked veery well in this instance. What we don't want is... thoughtful answer ... OP oh that wasn't the question in my mind. In this case that didn't happen

The Retrieve/ID mapping tool from Uniprot allows you to convert from Uniprot to entrez ID and vice-versa. It is available here. Just select the appropriate values in the 2.Selection options panel.

With thanks from me as well.

Thank you very much!


text:
How to reduce the occupied RAM when you are dealing with a very sparse matrix in a single-cell Experiment in R?
I'm dealing with a very large and sparse dataset and the first issues I met occurred when I tried to use quickCluster that reported me this error:
                'cannot allocate vector of size 156.6 Mb'

So, given that I cannot wait to change the RAM of my computer and I can't afford to use a cluster, I want to rely on some other strategies like some package that would allow me to handle sparse matrices. I'm thinking about sparseM but given that I don't know well this package I'd like to know how to shrink the ram allocation for these kind of matrices. Any suggestion will be very appreciated!
Link to the dataset

This is with the scran package, yes? Can you update your tags to mention this? To my knowledge scran doesn't support sparse matrix input, which will be a bit of a show-stopper unless you can get access to a system with more RAM (maybe you can pay a few bucks for an AWS instance?).

Yes it is. So you are telling me to utilize AWS?This could be a problem since it is a thesis and I do not want to pay.

You could go to GCP which is cheaper than AWS (Google Cloud). How many cores do you have ?

4 cores I presume

Ah, looks like I can't even procrastinate on StackExchange anymore without seeing work-related stuff. Oh well.
Anyway, the other answers and comments are way off. scran has supported sparse matrices for years, ever since we switched over to the SingleCellExperiment class as our basic data structure. quickCluster does no coercion to dense format unless you tell it so explicitly, e.g., with use.ranks=TRUE (in which case you're asking for ranks, so there's little choice but to collapse to a dense matrix).
You don't provide an MWE or your session information, but this is how it rolls for me:
# Using the raw counts in the linked dataset. Despite being
# called a CSV, it's actually space delimited... typical.
library(scater)
mat <- readSparseCounts("GBM_raw_gene_counts.csv", sep=" ")

# Making an SCE just for fun. Not strictly necessary for
# this example, but you'll find it useful later.
sce <- SingleCellExperiment(list(counts=mat))

library(scran)
system.time(clust <- quickCluster(sce))
##   user  system elapsed
##  3.170   0.174   3.411

This is running on my laptop - 16 GB RAM but I'm definitely not using all of it. I only go full throttle when I'm working on some real data, e.g., the 300k HCA bone marrow dataset. Check out the book for more details.
Session info below, I don't know quite enough SO-fu to know collapse it.
R version 4.0.0 Patched (2020-04-27 r78316)
Platform: x86_64-apple-darwin17.7.0 (64-bit)
Running under: macOS High Sierra 10.13.6

Matrix products: default
BLAS:   /Users/luna/Software/R/R-4-0-branch/lib/libRblas.dylib
LAPACK: /Users/luna/Software/R/R-4-0-branch/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets
[8] methods   base

other attached packages:
 [1] scran_1.16.0                scater_1.16.1
 [3] ggplot2_3.3.2               SingleCellExperiment_1.10.1
 [5] SummarizedExperiment_1.18.1 DelayedArray_0.14.0
 [7] matrixStats_0.56.0          Biobase_2.48.0
 [9] GenomicRanges_1.40.0        GenomeInfoDb_1.24.2
[11] IRanges_2.22.2              S4Vectors_0.26.1
[13] BiocGenerics_0.34.0

loaded via a namespace (and not attached):
 [1] beeswarm_0.2.3            statmod_1.4.34
 [3] tidyselect_1.1.0          locfit_1.5-9.4
 [5] purrr_0.3.4               BiocSingular_1.4.0
 [7] lattice_0.20-41           colorspace_1.4-1
 [9] vctrs_0.3.1               generics_0.0.2
[11] viridisLite_0.3.0         rlang_0.4.6
[13] pillar_1.4.4              glue_1.4.1
[15] withr_2.2.0               BiocParallel_1.22.0
[17] dqrng_0.2.1               GenomeInfoDbData_1.2.3
[19] lifecycle_0.2.0           zlibbioc_1.34.0
[21] munsell_0.5.0             gtable_0.3.0
[23] rsvd_1.0.3                vipor_0.4.5
[25] irlba_2.3.3               BiocNeighbors_1.6.0
[27] Rcpp_1.0.4.6              edgeR_3.30.3
[29] scales_1.1.1              limma_3.44.3
[31] XVector_0.28.0            gridExtra_2.3
[33] dplyr_1.0.0               grid_4.0.0
[35] tools_4.0.0               bitops_1.0-6
[37] magrittr_1.5              RCurl_1.98-1.2
[39] tibble_3.0.1              crayon_1.3.4
[41] pkgconfig_2.0.3           ellipsis_0.3.1
[43] Matrix_1.2-18             DelayedMatrixStats_1.10.0
[45] ggbeeswarm_0.6.0          viridis_0.5.1
[47] R6_2.4.1                  igraph_1.2.5
[49] compiler_4.0.0

Not a direct solution but some workarounds:

As far as I know, Seurat can work with sparse matrices.

The particular function of scran that you are using eats up quite some memory. I believe it is needed for the "normalization" step (that is how I used it anyway). Whereas the scaling normalization performed by this function is superior to crude "log normalization", you can give a try with the latter, which is far less computationally-intensive (does not do clustering). Seurat, once again, can help with this.

You can downsample your data to the extent that it fits your into your RAM.

You can give a try with Python. More and more single cell packages are written in Python, to some extent because of the problem you have experienced. For example Scanpy output is comparable to that of Seurat, although I am not sure if you can use scaling normalization with Scanpy.

I solved the issue. Anyway I will try what you suggests also because it is instructive. Thank you very much!

Please let the community know how have you solved the issue, this is important for people who will encounter the same problem. You can just "answer" your own question and "accept" the answer.

Essentially you've hit a RAM bottleneck and the calculation will slow to zero, or in this instance refuse to go forward. The way to do this normally is to parallelize the calculation across the cores of your machine. This will likely remove the RAM bottlenack, dont ask me the computer archetectural reasons for why it works - but its works.
However, my knowledge of R is minimal. I wouldn't know how to parallelize an R calculation. It is certainly doable in Perl and Python, but the calculation needs to be written to ensure parallelisation.
The other way in is to reconfigure your calculation to remove sparse matrices OR find someone doing NGS where they've confgured their machine around heavy RAM.
Looking at your calculation I don't quite understand why you need to use a specific package, it looks like unsupervised machine learning, PCA - tSNE that sort of thing and you don't need a given package to do that you just need to vectorise the inputs. If you worked out the statistical components of Scran then there are a few extremely strong R statisticians/bioinformaticians on site who would have zero problem in replicating this within a few lines of code. Its not hard in Python's Sci-kit learn either. At a guess they perform PCA and resolve it via tSNE and this give nice clear clusters.
GCP is free for 3 months so it will cost you zero in context to a single calculation.

Even PCA gives problem if it ingests this sce-object. Anyway I solved changing the computer with one more performative. Very good to know that GCP is free for three months. Thank you!

I prefer @haci 's work around, thanks are not a cool way to say thanks.but anyway

good suggestions @Michael. It's not a big dataset, 3500+ cells, so you can do PCA or tSNE on the dimension reduced data. I don't see the OP's calculation though, or what is the larger question..

.. me neither but they are doing Plasmodium (malaria, P. falciparum?) so there will be 24-hour periodicity in the RNA expression depending on what in vitro assay they are doing. That would make it a bit complicated. They also appear to do martial arts - its a tough field malaria ;-)

Not a direct solution to your problem because it does not use SCE. But if you (and others who ended on this post) are facing the RAM bottleneck when analysing single-cell data, then Scarf might be useful.
Here is the Github link:
https://github.com/parashardhapola/scarf
Documentation:
https://scarf.readthedocs.io/
Paper:
https://www.nature.com/articles/s41467-022-32097-3


text:
Comparing gene abundances between metagenomes
My workflow until now:
Find fragments of a marker gene in unassembled metagenomes > download and assemble metagenomes > recover the gene neighborhood / gene set of interest
Right now I have a rough estimate of how abundant these genes are by the 'depth' noted on the assembled contigs (I use MEGAHIT for assembly). I was wondering if there's a more thorough/proper way to do this. I would like to compare the abundance of specific genes between a) samples in the same study, and b) different studies. I imagine that the size of the individual metagenomes should be considered in both cases, but point b) might add additional difficulties such as different sequencing techniques. I would appreciate your insights.

I would avoid using assemblies to answer this question, as there's no guarantee that you will be able to assemble your genes of interest; you can however estimate their abundance even if they are relatively rare.
How I understand your question as being one of estimating the abundance of either some specific genes (e.g. butyrate metabolism genes) or all genes in a microbial community across multiple samples for comparative purposes. In other words, not 16S or marker gene analysis for the purposes of estimating organismal abundance, which is a rather different problem (though in that case I would still not use an assembly).
A more standard workflow is:

align metagenomic reads against some existing database of genes annotated appropriately.
estimate the number of reads aligning against some gene or ortholog in some way (using for example KEGG Orthology or similar).
use counts from (2) as input to some statistical procedure, possibly summarizing across functional categories.

Some examples of how this has been done are here, here, here. I am sure that there are more recent/relevant references but I haven't been following the field closely in the last few years.


text:
Comparing aligned amino acids to codon
I have a set of 4 amino acids which i have aligned,but i want to compare then with respective  nucleotide sequence in terms of change codon level.
I have the alignment file where i want to take only for the C-terminal domain region region only the CTD starts with amino acids “NITNLC” till it ends with “HAPATV” and rest of the sequence i dont want to take.
In short i want to take the part sequence which starts from "NITNLC" to "HAPATV" then i would like to compare what are the changes in the nucleotide sequence at codon level.

To start with I have used the seqinr library to generate nucleotide alignment from aligned protein sequence.
reverse.align {seqinr}

Then the next step that i want to compare each of my amino acid sequence to my reference sequence which in my case is P0DTC2to rest of the sequence for example I want to compare P0DTC2 to K9N5Q8to the above mentioned CTD region which starts from NITNLC and ends with HAPATVand find what are the codon level changes if there is any I would like to report both in amino acid level and codon level.

First part I guess my approach is correct. Second part I'm not sure how to proceed i guess it's much more than simple parsing I guess!!.
Any help or suggestion would be really appreciated ,and if there exist a R based solution that would be really welcome.
Files
Amino acid file nucleotide file Aligned amino acid file Reverse alignment nucleotide file using seqinr

To set the ball rolling, I tried to implement it manually with your files. Overall, there are three steps:

Locate NITNLC (or HAPATV) in protein P0DTC2.
Locate mismatched amino acids between proteins P0DTC2 and K9N5Q8 within the range from step 1.
Print the amino acid and the DNA codon of both proteins from step 2.

It works, but only for the first 60 amino acids. I wonder if it's due to the amino sequences in AMINOOO_seq_removed.fasta repeat itself every 60 acids. But why?
#For example, the first three lines of protein P0DTC2
>P0DTC2
MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHSTQDLFLPFFS MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHSTQDLFLPFFS
NVTWFHAIHVSGTNGTKRFDNPVLPFNDGVYFASTEKSNIIRGWIFGTTLDSKTQSLLIV NVTWFHAIHVSGTNGTKRFDNPVLPFNDGVYFASTEKSNIIRGWIFGTTLDSKTQSLLIV
NNATNVVIKVCEFQFCNDPFLGVYYHKNNKSWMESEFRVYSSANNCTFEYVSQPFLMDLE NNATNVVIKVCEFQFCNDPFLGVYYHKNNKSWMESEFRVYSSANNCTFEYVSQPFLMDLE
...

Step 0: Read the files.
library('seqinr')

#align file containing protein sequences
count_added <- read.alignment('count_added_.clustal_num', format='clustal')
names(count_added$seq) <- count_added$nam

#DNA sequences
rev3.aln <- read.alignment('rev3.aln', format='fasta')
names(rev3.aln$seq) <- rev3.aln$nam

Step 1: Locate NITNLC (or HAPATV) in protein P0DTC2. Again, there are two repeats of NITNLC 60 acids apart (at 807 and 870).
locate <- function(seq, find)
  {address <- gregexpr(paste(strsplit(find, '')[[1]], collapse='[^a-z]*'), seq)
  #substr(seq, address[[1]][1], address[[1]][1]+attr(address[[1]], 'match.length')[1]-1)

  return(list(start=as.numeric(address[[1]]),
             end=as.numeric(address[[1]] + attr(address[[1]], 'match.length') - 1)))
  }

locate(seq=count_added$seq[['P0DTC2']], find='nitnlc')
#start
#807 870
#end
#812 875
locate(seq=count_added$seq[['P0DTC2']], find='hapatv')
#start
#1212 1279
#end
#1217 1285

Step 2: Locate mismatched amino acids between proteins P0DTC2 and K9N5Q8 within the range from step 1. Instead of using the range 812 - 1279, I chose range 1 - 20 for demonstration purposes.
compare <- function(seq1, seq2, after=0, before=100000)
  {seq1_ = strsplit(seq1, '')[[1]]
  seq2_ = strsplit(seq2, '')[[1]]
  ind = which(seq1_ != seq2_ & grepl('[a-z]',seq1_) & grepl('[a-z]',seq2_))
  ind = ind[ind>after & ind<before]
  #seq1_[ind[1]]
  #seq2_[ind[1]]

  return(ind)
  }
compare(seq1=count_added$seq[['P0DTC2']], seq2=count_added$seq[['K9N5Q8']], after=1, before=20)
# [1] 5  7  8  9 10 13 14

#protein comparison
#K9N5Q8      MIHSVFLLMFLLTPTESYVD
#P0DTC2      ----MFVFLVLLPL------
#<mismatch>      5 7890  34

Step 3: Print the amino acid and the DNA codon of both proteins from step 2. Note that ind is based on the align file.
print_amino_codon <- function(ind, seq, seq_gene)
  {locate_amino <- gregexpr('[a-z]', seq)[[1]]
  if (!ind %in% locate_amino) return(NA)
  ind2 = match(ind, locate_amino)

  return(c(amino=substr(seq, ind, ind), codon=substr(seq_gene, ind2*3-2, ind2*3)))
  }

codon(ind=5, seq=count_added$seq[['K9N5Q8']], seq_gene=rev3.aln$seq[['K9N5Q8']])
#amino codon
#  "v" "atg"
codon(ind=6, vseq=count_added$seq[['K9N5Q8']], seq_gene=rev3.aln$seq[['K9N5Q8']])
#amino codon
#  "f" "ttc"
codon(ind=7, seq=count_added$seq[['K9N5Q8']], seq_gene=rev3.aln$seq[['K9N5Q8']])
#amino codon
#  "l" "ttg"

codon(ind=5, seq=count_added$seq[['P0DTC2']], seq_gene=rev3.aln$seq[['P0DTC2']])
#amino codon
#  "m" "atg"
codon(ind=6, seq=count_added$seq[['P0DTC2']], seq_gene=rev3.aln$seq[['P0DTC2']])
#amino codon
#  "f" "ttt"
codon(ind=7, seq=count_added$seq[['P0DTC2']], seq_gene=rev3.aln$seq[['P0DTC2']])
#amino codon
#  "v" "gtt"

#protein comparison
#K9N5Q8      MIHSVFLLMFLLTPTESYVD
#P0DTC2      ----MFVFLVLLPL------
#<mismatch>      5 7890  34
#<print>         ^^^

#K9N5Q8 gene codon
#gtg ttt cta ctg atg ttc ttg tta aca
#                ^^^ ^^^ ^^^
#P0DTC2 gene codon
#atg ttt gtt ttt ctt
#^^^ ^^^ ^^^

Perfect @RyanSYKwan its alot of code

thanks a lot this will be really really helpful for my future work

Glad it helps! I learn a lot from you too.


text:
Python: How to write duplicate sequences removed from fasta file to new file
I currently am using this code to remove duplicate sequences from the fasta file. However, I would also like to write a new file with only the removed duplicates as well as a count for how many times they appear. Please let me know if you have any suggestions, I am hoping to keep it written in one py file, but I am not opposed to running a second file to do so. I attended to write the appended strings to a new file but it looked slightly off...
from Bio import SeqIO
import time

start = time.time()

seen = []
records = []

 for record in SeqIO.parse("b4r2.fasta", "fasta"):
   if str(record.seq) not in seen:
       seen.append(str(record.seq))
   records.append(record)

#writing to a fasta file
SeqIO.write(records, "no_dupes_b4r2.fasta", "fasta")
end = time.time()

print(f"Run time is {(end- start)/60}")
```

Hi @User this has been nicely answered from different perspectives. Could you please send off one of the answers to this question as "accepted". The site has certainly helped here, so to reciprocate that help, this certainly helps.

Edit: user @RamRS has a more robust answer for actually using this. Please use that. Leaving answer up as OP appears to be interested in implementing as a learning activity.
Please let us know why this is not doing what you intend.
I haven't checked it, but something like this would work:
from Bio import SeqIO
import time
from collections import defaultdict

start = time.time()

seq_counts = defaultdict(0)
records = []

 for record in SeqIO.parse("b4r2.fasta", "fasta"):
   seq = str(record.seq)
   if seq in seq_counts:  # python is smart about dict keys
       records.append(record)

   seq_counts[seq] += 1

#writing to a fasta file
SeqIO.write(records, "no_dupes_b4r2.fasta", "fasta")
end = time.time()

# print the counts
for seq in seq_counts:
    print(seq_counts[seq], seq)

print(f"Run time is {(end- start)/60}")

It will be a bit slow for large fasta files and I am sure that some people will be very upset that I'm using a defaultdict for this.
The code that you have would probably not work as you would want the records.append() call to be only run in the if statement testing for non-membership.
There is absolutely no reason to have to write more than one script for this.

There is absolutely no need to write even this script when tools like seqkit rmdup already exist.

Hi @RamRS any chance you could put this as an answer? It would be helpful for my own personal records if nobody else.

Thank you! I've seen some other scripts that can do similar things but I wanted to make my own... :) Appreciate it a lot.

i am indeed upset that you use collections.defaultdict, would prefer collections.Counter ;) e.g. seq_counts = Counter(record.seq for record in SeqIO.parse(...

@RamRS I'll note that seqkit rmdup does not appear to fully satisfy the asker's question, which asked for a specific count of each sequence.

@Chris_Rands yeah totally, i just know the defaultdict syntax by heart and i'd have to go look up how to use a counter again...

@MaximilianPress -D and -d do that. If OP wants to save count against each sequence (as opposed to against each sequence identifier), the file would be unreliable as it would need to be a FASTA with an annotation added in the header. With the -D, OP should be able to use seqkit grep and seqkit duplicate to get their custom file.

@RamRS all right then, that is what they should do if they're using this in earnest. On the other hand, it does seem like they want to implement it themselves as an exercise, so I'll leave my answer up for now with an edit.

Agree with you 100%. This is great as a learning exercise (who among us did not write our own fasta deduplicator?), but these robust tools, as you rightly call them, are better tested so they catch edge cases that our code might not look out for.

Please do not reinvent the wheel - seqkit rmdup does what you want.
From the website:

rmdup
Usage
remove duplicated sequences by id/name/sequence

Usage:
  seqkit rmdup [flags]

Flags:
  -n, --by-name                by full name instead of just id
  -s, --by-seq                 by seq
  -D, --dup-num-file string    file to save number and list of duplicated seqs
  -d, --dup-seqs-file string   file to save duplicated seqs
  -h, --help                   help for rmdup
  -i, --ignore-case            ignore case

How does a -D file look like, you ask:
$ zcat hairpin.fa.gz \
    | seqkit rmdup -s -i -o clean.fa.gz -d duplicated.fa.gz -D duplicated.detail.txt

$ cat duplicated.detail.txt   # here is not the entire list
3   hsa-mir-424, mml-mir-424, ppy-mir-424
3   hsa-mir-342, mml-mir-342, ppy-mir-342
2   ngi-mir-932, nlo-mir-932
2   ssc-mir-9784-1, ssc-mir-9784-2

Hi @RamRS perfect thank you very much, that is really helpful


text:
Coding vs non-coding DNA length
I am trying to calculate the total exonic length (in bp) in order to see where the coding-noncoding ratio of roughly 1% to 99% comes from.
We know all chromosomes total about 3 billion base pairs. So Is it correct to assume that coding regions would total about 30 million base pairs? I am looking at the numbers in the sources below, but nothing there adds up...
Take numbers from NCBI annothation for example; specifically

no matter how i count i can't get the expected ratio... What am  doing wrong?
Thank you!
https://www.encodeproject.org/report/?type=Gene&organism.scientific_name=Homo%20sapiens&searchTerm=gene&field=symbol&field=notes
https://archive.vn/20130414235101/http://useast.ensembl.org/Homo_sapiens/Location/Chromosome?r=1

cross-posted: https://www.biostars.org/p/447675/

There are a few different answers provided here: https://biology.stackexchange.com/questions/48110/how-is-the-size-of-a-gene-defined

Thank you, but this is a different question; my aim is to get the said ratio given the feature numbers, in my case the length of coding vs non-coding regions. The original response helped me by pointing out that i should focus on exons, rather than the entire gene region. So i searched and found some definitive numbers for exon length, and calculated in different combinations, but they still aren't adding up to the 99%-1% ratio. So i'm wondering if i'm missing something else here. Or whether i'm looking at the wrong features entirely.

They are not adding up because you need exclude non-coding exons (non-CDS) and exclude overlapping features. The discussion I posted shows how to get coding regions. The only difference is that they are trying to get the average length and you need total length.


text:
Question about Using MDAPackmol (packmol in python)
I would like to use the python wrapper for packmol (MDAPackmol) for my research. In order to do so, I have tried to work with the example code and pdb files provided on GitHub so I can gain an understanding of the how to use the wrapper. However, when I run the code, I receive the following error message: ValueError: invalid literal for int() with base 10: '' .
Here is the complete trace back.

The error originates at start.resname being a zero-length string. It seems to be a problem with the resname having a number and a letter and thus is a problem with Universe (the object returned by load_packmol_output). Thus, there is a problem with the residues. Any advice or insight would be greatly appreciated. Thanks!

H @bioinf99 thanks for your post. The error is simple the parser is expecting an integer and not getting one. Tracing it on the information given, is not going to be easy. The location doesn't help much. It ain't packmol thats the problem. Maybe someone here has used the wrapper ... It could lots of things ...I ain't optimistic here, but lets see. My best guess is its the 'argparse' command and you're not entering what it wants to see. In that case you'd need to look at the opening code. If that doesn't mean much to you I'd contact the developers. Just my views however.

Thanks for your reply. I agree with your assessment. I have looked at the source code but am unsure about making changes directly to it. Unfortunately, I have emailed the developer, but he has not answered.

You'd need to send it through an IDE with breakpoints (or at least that's how I'd do it, to look at all the variables and their parsing to the package) and thats lots of work. Someone on here might know, or the developer might be on site. Consider Snakemake if this doesn't work out?? Snakemake is well supported (I don't use it BTW)..

Thanks for the ideas. Definitely will look into it if there isn't an easier solution.

One thing worth checking is if test files provided crash similarly. Also, about the email radio silence, a package with 7 stars could be one that is not really intended for the general public. An alternative is raising issues on GitHub —I for one answer those quicker as watchers get notified too and is public.

Thank you for your comment. The test file does not appear to produce any errors. Thank you for your insight regarding raising the issues on GitHub. My lab really wants this to work.

Sorry, I just re-read your Q and I think you've inadvertently misdirected everyone. The first arrow points at the function/class call not at a specific argument. The error is from start.resname being a zero-length string —not sure why resname is a letter + number. It is a property of the Universe object returned by load_packmol_output. So the packmol call seems to have worked, but the residues have some oddity. So you might want to alter your Q...

Thanks. So don't worry about the instructions? That would make sense, because it is what was shown in the documentation. It's actually a problem with MDAnalysis? Thanks again.

Sorry for another question, but I am new to working through these types of errors. I was wondering how you knew it started at start.resname. Thanks!


text:
RNA seq .counts.txt to bigwig conversion
I've downloaded neuronal RNA seq data from GEO. The files are in .counts.txt format. Could I convert them to bigwig? If so, how?

Hi @Samir ... you've got the files from here https://www.ncbi.nlm.nih.gov/geo/ cool, could you just fill in details about viewing and bigWig format for the rest of us?

Hi Michael. I want to load RNA seq data into HiGlass, a browser that I use to look at my HiC data. I have previously loaded RNA seq bigWig files into HiGlass. I just don't know if there is a way to load the text file format into HiGlass, so I wanted to know if I can convert the format to bigWig.

can you provide the GEO accession ? Sometimes they do provide bigwig but i suspect it's not very useful for RNAseq unless there's something specific u need

Hi, this is the GEO accession: GSE102838

No, count-level data doesn't have the base-level resolution that is required in order to be usefully visualised on a genome browser. You need a BAM file (or raw FASTQ files that can be used to generate a BAM file) in order to do that.
"No" isn't completely correct: it'd be possible to take gene counts and information about gene annotation boundaries to shoehorn them into a BedGraph / BigWig file, but whatever comes out the other end of that isn't going to be a good representation of reality.


text:
Trying to retrieve sra data using sra explorer
I am trying to use the SRA explorer tool to download data from ENA. Whatever geo accession number I give it however, it will return no results. It only seems to return results for the suggestions it already has. I tried for example searching for GSE111797 by typing it in the search tool but didn't get anything.

It would be helpful to describe the question in more detail, like what you did more specifically. Also, since you are using a third-party tool, it might be easier to reach out to the developers themselves.

Try the corresponding SRP accession (SRP135519). It seems to work for me now.

@vkkodali yup that did it


text:
Error pasing the following line in pdb
I converted an .sdf file containing around 50,000 anti-viral agents to individual .pdb files using OpenBabel. However, when I tried opening the .pdb files with Autodock I kept getting the error below. What do you think is causing the error.
Thank you for your time and patience.
Error Message:

“ Error pasing the following line in pdb:
HETATM 1 C UNL 1 52355.00975659.051 0.000 1.00 0.000 ."

You have a firewall issue, it is needing to connect outside your computer.

I have added autodock.exe and autogrid.exe as exceptions to my firewall and yet I am still facing the same issue.
Here’s is exact the error message,
“ Error pasing the following line in pdb:
HETATM 1 C UNL 1 52355.00975659.051 0.000 1.00 0.000 ."
What do you think I should do to overcome the error. Thank you.

Well at least it is now an incomprehensible error - but it still doesn't work :-(  .. some of the pdb gurus are needed here, its a very ancient format BTW.

A PDB file has strict spacing rules as documented here.
The coordinates should be:
31 - 38       Real(8.3)      x             Orthogonal coordinates for X.
39 - 46       Real(8.3)      y             Orthogonal coordinates for Y.
47 - 54       Real(8.3)      z             Orthogonal coordinates for Z.

So trying this on your line:
'HETATM 1 C UNL 1 52355.00975659.051 0.000 1.00 0.000'[30:39]

I get 9.051 0., which is not a number.
The whole line is formatted incorrectly however.
So as you say you are using Autodock why not pass a mol2 to prepare_ligand4.py instead which is an acceptable format?
Alternatively, you could convert your multi-entry sdf to pdb files via RDKit in Python:
suppl = Chem.SDMolSupplier(sdfile)
for mol in suppl:
      # mol = Chem.AddHs(mol) # has hydrogens?
      # AllChem.EmbedMolecule(mol) # is 3D?
      AllChem.ComputeGasteigerCharges(m) # sdfs do not have partial charges
      Chem.MolToPDBFile(m, f"{mol.GetProp('_Name')}.pdb")

Hello Mr. Ferla,
Thank you for getting back to me.
Is there any option in Open Babel that will let me make the required changes so that the spacing rules aren't violated.
Best Regards,
Shahbaaz.

@Shahbaaz , amended.
Parenthetically no need for formalities. Also @ plus username is generally the preferred way to address folk here, you might get titles  wrong ;) .

Also, this other Q about Autodock protocol may be useful: https://bioinformatics.stackexchange.com/questions/11471/protein-ligand-docking-how-to-convert-protein-pdb-to-protein-maps-fld/11479#11479

Thanks a ton for the help

Good morning @Matteo Ferla
Can you briefly summarize how I can use prepare_ligand4.py or RDKit. I have looked for tutorials online but couldn't come across any. Thank you.

First, prepare_ligand4.py is for Autodock 4, not VINA. The two are different and fill different roles: complicated but adaptable vs. easy but blackbox. I am assuming you are using the former as you did not specify VINA. Here is again the link to documentation.

I need help with the following line,
"pythonsh [yourpath or .]/prepare_ligand4.py -l ligand.pdb [options]"
considering the path for mg tools is the following, "C:\Program Files (x86)\MGLTools-1.5.6" how would I alter the above line including the mg tools file path in it

If you use finder (whatever it is called on a Windows) and explore the stuff within C:\Program Files (x86)\MGLTools-1.5.6 you'll find the script in somewhere like autodock_tools/utilities

I have tried everything to split the multi entry sdf file to pdb files but have failed. Could you split the sdf file to individual files for me as I have been unable to do so despite multiple attempts at it.
Thanking you in anticipation,
Dr. Hussain.

If it is just for training, why not redock ligands from already bound structures one at the time. There are lots of inhibitors to the COVID protease in the PDB and the ligand can be downloaded as sdf, pdb or mol2. Make an apo, place the ligand (don't keep the original there) and then you can see what RSMD you get versus the actual solved one.
Once there you can try your hand at virtual library screening.


text:
How can I get unmatched reads for defective genomes analysis using bwa and samtools?
I am trying to follow the Beauclair et al. paper (free version here https://rnajournal.cshlp.org/content/24/10/1285.long) for identifying defective genomes using their DI-tector program.
According to figure 2
Before using the program itself, it is highly recommended to filter my raw read using bwa and then samtools
Unfortunately I have no idea where to start - I have browsed official bwa and samtools sites, I have installed them using the following commands:
sudo apt-get update -y
sudo apt-get install -y bwa

I have indexed my reference genome with the following command (as I read it is required before the analysis)
bwa index sequence.fasta

Then, I have installed samtools
cd ..
wget https://github.com/samtools/samtools/releases/download/1.9/samtools-1.9.tar.bz2
tar -vxjf samtools-1.9.tar.bz2
cd samtools-1.9
make
sudo apt install samtools

I wonder what will be the further steps to get these "UNMATCHED READS at the end of step i", since Di-tector program takes virus_reference genome (I have it) and these unmatched reads.
First of all I want to align my .fastq files vs GCF_000001405.39_GRCh38.p13_genomic.fna - what will be the output file? If I got it right - I need this output file to filter with samtools, am I correct?
If I am not mistaken, the output after samtools will be the .sam file that I will have to put into a DI-tector - correct?

If you want to keep track of the unmapped reads, I recommend you to use Bowtie2 instead of bwa.
Start by building an index:
bowtie2-build genome.fa genome.fa

Then map your reads, using arguments to save unmapped reads in a specific file like this:
bowtie2 -x genome.fa -1 reads_1.fq -2 reads_2.fq -S bowtie2_mapping.sam --un-conc-gz unmappedreads_%.fq.gz

It will write the alignment as well as unmappedreads_1.fq.gz and unmappedreads_2.fq.gz that you may use in DI-tector.

My thoughts exactly. Personally I don't know why a paper in 2018 would want to use BWA over Bowtie2. Any ideas?

after putting a glance at the article, it is no really clear how they switch from bwa mem to bwa samse to align reads. Maybe it has something to see with the flags in the output .sam file, required for DI-tector proper running. But it is only a supposition...

Thank you. But what about samtools step? I mean, my read has to be aligned with bwa (or bowtie2 - great I will try that) FIRSTLY with the human genome, then SECONDLY with virus genome and only after these two steps it has to be filtered through SAMTOOLS and only after that I have to have the unmapped reads that I can put into a DI-tector.

I think it depends on the format of input data required for DI-tector. Bowtie2 will return .fastq files. In the article, they used samtools with -Bs -f 4 options probably to return unaligned reads into a .sam with specific CIGAR labels.

@thomas duge de bernonville OK, I did bowtie2-build GCF_000001405.39_GRCh38.p13_genomic.fna human_genome.fa first, and the output is 4 files human_genome.fa.1.bt2, human_genome.fa.2.bt2, human_genome.fa.3.bt2, human_genome.fa.4.bt2. I wonder how can I map my reads using these 4 files while in the command you have proposed you have only genome.fa Anyways I have already tried this bowtie2 -x human_genome.fa /mnt/e/nastya/SLX066-02/B-dVMV-RIG-1/B-dVMV-RIG-1_ACAGTG_L008_R1_001.rc.fastq.gz -S bowtie2_mapping.sam --un-conc-gz unmappedreads_%.fq.gz and it gave me an (ERR) bowtie2 exited

try with bowtie2 -x human_genome.fa -U /mnt/e/nastya/SLX066-02/B-dVMV-RIG-1/B-dVMV-RIG-1_ACAGTG_L008_R1_001.rc.fastq.gz -S bowtie2_mapping.sam --un-conc-gz unmappedreads_%.fq.gz, adding the -U argument to indicate the single end read file

18035544 reads; of these:   18035544 (100.00%) were unpaired; of these:     17979950 (99.69%) aligned 0 times     55594 (0.31%) aligned exactly 1 time     0 (0.00%) aligned >1 times 0.31% overall alignment rate this is what I have now. I also have unmappedreads_1.fq and unmappedreads_2.fq archives of 1Kb each.. Is that ok? Why are there 2 files? Which one can I use with DI-tector? I also got a bowtie2_mapping.sam file which is about 3Gb. Thanks!

This is because you asked for it --un-conc-gz unmappedreads_%.fq.gz. But I don't know if your input data is indeed paired. What is the size of each unmapped*fq?  What input is required for DI-tector?

Thanks again, and sorry for taking your time. There are two unmapped reads and the size is 1Kb each. In fact, they are empty. DI-tector requires unmatched reads in .fastq format + Virus_ref genome (I believe it's an indexed .fasta viral genome)

ok I think it is my fault, you do not have paired-end reads but single end. So the correct option should be --un-gz, instead of --un-conc-gz, as described in the bowtie2 manual.


text:
How to aggregate multiple MAF files together?
I have cancer genomic data (tumor/normal whole exome sequencing) from 50
patients that received the same type of treatment, half of whom responded. These come in the form of 50 .maf files, along with a supplemental file that, along with some other fields, has the Response (Responder vs Non-Responder) field. My question is how to aggregate all of this together so that I can perform some statistical tests on the data. I have a clever way (taken from this blog post) of reading in the 50 .maf files and combining them, but I wonder if this is an appropriate approach.
sample_info <- readr::read_tsv(file = "path/to/sample_info/sample-information.tsv")
maf_files <- fs::dir_ls("path/to/mafs/")
patients_data <- maf_files %>%
  purrr::map_dfr(read_tsv, col_types = list(Chromosome = col_character()))

My thought then was to just dplyr::left_join() the sample_info with patients_data like
patients_data_final <- patients_data %>% dplyr::left_join(sample_info, by = c("Tumor_Sample_Barcode", "Matched_Norm_Sample_Barcode"))

For clarity, here are the column names of both dataframes
> colnames(patients_data)
 [1] "Hugo_Symbol"                 "Chromosome"                  "Start_position"
 [4] "End_position"                "Variant_Classification"      "Variant_Type"
 [7] "Reference_Allele"            "Tumor_Seq_Allele1"           "Tumor_Seq_Allele2"
[10] "Tumor_Sample_Barcode"        "Matched_Norm_Sample_Barcode" "Protein_Change"
[13] "t_alt_count"                 "t_ref_count"

> colnames(sample_info)
[1] "Patient_ID"                     "Tumor_Sample_Barcode"           "Matched_Norm_Sample_Barcode"
[4] "Response"                       "Silent_mutations_per_Mb"        "Nonsynonymous_mutations_per_Mb"
[7] "Mutations_per_Mb"

My task is to find genes that "are enriched in patients who either
responded or not." So as a supplemental question, if anyone has suggestions on which statistical test to use (or how to go about deciding), I'd appreciate that as well.
PS: I am aware of the maftools package which probably has an easy solution to this, but unfortunately my machine is old (Late 2011 MacBook Pro) and unable to run it. (Old Mac --> Can't update OS --> Cant update version of R --> Can't install necessary packages to run maftools)

Ask someone with a newer computer for help, or fire up an AWS instance with Ubuntu and do it there. You will face a lot of difficulties in bioinformatics if your machine is >5 years old and you don't have access to any HPC. Old notebooks can work well as terminals but not as workhorses.


text:
Encoding a 15bp DNA sequence into a shorter than 15 digits number
I have a series of 15bp DNA sequences from single-cell barcodes and I would like to encode them into a shorter than 15 digits number so that each unique 15bp DNA sequence will create a unique numeric value. What's the best way of achieving this?

What you are asking is called a "hash" (wiki).
In python, get your oligo as a string and use the standard library method hash(oligo), this is a number that gets used for a lot of operations. It is however always a 18+1 digit number in decimal, so a different solution is needed.
There are four nucleobases, so binary conversion from scratch is easy and best option and each solution will be unique and small. Each base takes up 2 bits and you concatenate these and convert them to hex or decimal —the length should be 10 decimal digits. One issue is when you have a oligo of a different length which start with adenosine or not.
One issue with this is the order. I am not sure where it comes from but the bases are nearly often expressed in the order ATGC. So...
mapping = {'A': '00', 'T': '01', 'G': '10', 'C': '11'}

It is satisfying how the first digit is the W/S, while the second is R/Y.
Now simply covert and add them compressed = functools.reduce(lambda acc,base : acc+mapping[base], oligo, '') and covert the string to decimal like so int(compressed, 2).

The problem is hash collisions, there is no guarantee that the compression is unique

Hash collision with the fake-hash would happen with two oligos that are the same but one has one or more extra adenosines. The oligo was said to be 15bp, so this does not apply. A workaround is to add a 1 at the start...


text:
Is there an "official" standard usage of Red and Green on differential gene expression heatmaps?
I have recently found myself making multiple heatmaps for visualization of differential gene expression results between replicates in two experimental conditions.
The default settings in the plotting software I use color:

Upregulated genes --> red
Downregulated genes --> green
Unchanged genes --> black

My intuition screams at me that this color scheme is reversed, and that upregulated genes should be green and that downregulated genes should be red. I have asked several labmates about their intuition regarding this color scheme, and have found about a 50/50 split in those that intuitively agree with each different scheme.
In my searching on this topic, I have found that several programs plot heatmaps with the same default colors as the program I'm using (red=upregulated, green=downregulated). Additionally, I have been told in conversation that this is essentially an accepted standard.
If this is indeed the convention, I will begin trying to rewire my intuition, but I want to confirm that it is before I do so. However I have not yet been able to find any "official" guidelines from agencies, publishers, authorities in the field, etc., specifying the use of red/green in heatmaps.
I do know that there are many suggestions for alternative heatmap color schemes, but I'm interested in finding an official answer for the red/green scheme specifically.
Edit: Several comments have stated that there is no such standard as described above. Therefore I am modifying the following question somewhat to:
Are there "official" guidelines from agencies, publishers, authorities in the field, etc., specifying the use of red/green in heatmaps?
Or conversely are there any “official” guidelines stating that there is no such specification and/or that any such selection of colors is acceptable?

No, there are no official guidelines.

If it is an accepted standard, that seems hard to believe.

There is no accepted standard on this and the somewhat common red/green coloring is absolutely terrible (and backwards).

Red/green is generally not a great colour scheme since it is a common colour blindness

Personally, I try to avoid heatmaps as much as I can. Usually, I find them confusing and I find it hard to mentally convert colour intensities to quantities, no matter what the colour scheme is. Having said that, I spent the whole of today dealing with heatmap.2. BTW, if you pin me down I'm on the side of upregulated=red (don't ask me why). Sorry, not a helpful comment...

No standard here, and I avoid plain red and green 1) because it has some tendencies that green is considered good/safe/trustworthy while red is the opposite and 2) there are red/green-blind people.

As others already pointed out, there is no official heatmap color scheme. I am not sure there is an official color scheme for any other types of plots either.
However, heatmaps seem to be far more standardized. Traditionally, heatmaps have been red-black-green as you originally stated. Most heatmaps from the microarray era and early RNA-seq are like that. You can have a heatmap for any type of data, but the majority is for gene expression. A lot of tools and tutorials were written back then, so they have that color scheme. From about 2008 to 2012, most people switched to the red-white-blue (or red-yellow-blue) color scheme. Switching green to blue avoids the color-blindness issue. Switching black to white is better for the ink if you print your papers and generally looks better on white background.
I originally found red as high to be confusing because I was thinking in financial terms (if the market is red, it's down). Now I think of it as a water faucet (red is hot, blue is cold).
Edit: In the last year or so, the viridis palette started to become widely used. It's far from the standard, but could become more popular in the future.

I really like this analogy, the Red is hot, Blue is cold will probably help me a lot! I agree that the "financial" use of red and green is probably the source of my switched intuition.

Although there is no standard as DeepTools' Devon Ryan said, there are still some habits. Usually red means upregulated, as you can see in the default color scheme in MATLAB's clustergram function. But you also want to make sure color-blinded people can see the difference. So another color scheme might be better, such as red=high while blue=low.
Anyway, if you can easily flip the color, using red to represent high is expected.

In the absence of any yet-identified official recommendations, I will update this answer with resources as I find them:
Example usage in tutorials:

EMBL Training: Red=Up, Black=NC, Green=Down

Program Defaults:

PATRIC Webiste: Red=Up, Black=NC, Green=Down

Guides:

Bio Turing Blog: Dos and don’ts for a heatmap color scale
sebastianraschka: A short tutorial for decent heat maps in R

Tools (by year):

(-2020) Web, R, CLI: MORPHEUS: Versatile matrix visualization and analysis software
(2019) Web, R, CLI: NOJAH: NOt Just Another Heatmap for genome-wide cluster analysis
(2017) Web, Jupyter: Clustergrammer, a web-based heatmap visualization and analysis tool for high-dimensional biological data
(2017) Web, R: shinyheatmap: Ultra fast low memory heatmap web interface for big data genomics
(2016) R: Complex heatmaps reveal patterns and correlations in multidimensional genomic data
(2016) Web: Heatmapper: Web-Enabled Heat Mapping for All
(<2016) GUI, Java, R: GENE-E
(2015) R: Heatmap3: an improved heatmap package with more powerful and convenient features
(2011) GUI, CLI: Gitools: Analysis and Visualisation of Genomic Data Using Interactive Heat-Maps


text:
How to anonymise a bam file
I have a bam file of, for instance, RNA-Seq, which contains patient-identifiable data in the form of Single Nucleotide Polymorphisms (SNPs) throughout.
I would like a technique to take this aligned bam file and replace all the non-reference bases with the reference bases at the equivalent positions (hard-clipping any soft-clipped bases while doing so), in order to have an anonymised bam file without the patient-identifiable information.
What is the best way to do this? I want the resulting anonymised bam file to have the reads in the same positions, for downstream tools that care about the alignment positions.

Hi @KraZug normally patient confidentiality is maintained through anonminity of name and personal information, but not to anonymize the DNA by ad hoc editting. No ethics committee would request that. Personally, I don't agree here.

In the event of a publication, normally you can submit to ega or dbgap https://www.ncbi.nlm.nih.gov/gap/docs/submissionguide/#:~:text=To%20comply%20with%20HIPAA%2C%20personally,submitted%20in%20dbGaP%20data%20files. and they have a set of guidelines on how u de-identify samples

You can simply replace the sequence and quality scores with a single * in each field. The CIGAR information will still be there, so it will be apparent where there were InDels and possibly mismatches (assuming your aligning uses = and X CIGAR operations), but since the actual sequence of these variants cannot be known this is quite anonymized.


text:
Nanopore microbiome data analysis 16s rRNA data
The steps I have done so far:

used guppy to do base calling from fast5 to fastq. In each barcode file I have 1960 fastq files which I have merged.
trimmed the fastq files using NanoFilt

Now what I see from various resources is they do de-novo assembly using tools like canu followed by minimap.
So is it needed to do de-novo assembly?
Or I can simply use my fastq files and process them using tools like kraken for taxonomic classification?
Any help or suggestion would be really appreciated


text:
Problem with merge data while trying to convert gene names in R
This question has also been asked on Biostars and StackOverflow
I've been trying to code (in R) a way to convert gene accession numbers to gene names (from RNAseq data). I've looked at all the related questions and tried to modify my code such, but for some reason it's still not working. Here is my code, where charg is a character vector of the gene accession ID's of the data set resdata:
charg <- resdata$genes
head(charg)

library(biomaRt)

ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl")

theBM <- getBM(attributes='ensembl_gene_id','hgnc_symbol',
      filters = 'external_gene_name',
      values = charg,
      mart = ensembl)

resdata <- merge.data.frame(resdata, theBM, by.x="genes",by.y="ensembl_gene_id")

Here's some output (where I'm struggling):
> head(charg)
[1] "ENSG00000261150.2"  "ENSG00000164877.18" "ENSG00000120334.15"
[4] "ENSG00000100906.10" "ENSG00000182759.3"  "ENSG00000124145.6"

> dim(theBM)
[1] 0 1

> head(theBM)
[1] ensembl_gene_id
<0 rows> (or 0-length row.names)

> dim(resdata)
[1] 20381    11
> resdata <- merge.data.frame(resdata, theBM, by.x="genes",by.y="ensembl_gene_id")
> dim(resdata) #after merge
[1]  0 11 #isn't correct -- just row names! where'd my genes go?

Thank you.

This is the code to get a look-up table to convert between Ensembl ID and HGNC:
ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl")
theBM <- getBM(attributes=c('ensembl_gene_id','hgnc_symbol'),
               filters = c('ensembl_gene_id'),
               values = gsub("\\..*", "", charg),
               mart = ensembl)

What Devon was posting is correct but misses a c() around the attributes values.
For further help please provide the content of resdata which you should always do when posting a question, since we cannot read minds. Does not work by the way is not a proper error description.
Once you have the output do:
resdata$genes <- gsub("\\..*", "", resdata$genes)

merge(x = theBM,
      by.x = "ensembl_gene_id",
      y = resdata,
      by.y = "genes")

Note that I had to go to that SE crosspost to get the content of resdata, this is not how this goes. Please post all relevant data up front in the future otherwise your questions might get downvoted and closed. Please also avoid cross-posting. if you provide proper information you usually get a good answer in time.
Edit: Just realized you also cross-posted this to Biostars even twice. Please stop this. I closed the Biostars posts and gave my two cents on this behaviour over there.

Truly, I apologize for any inconvenience my multiple posts caused, and it certainly will not happen again. Thank you an incredible amount for helping me, this worked beautifully :). Have a wonderful day.

Glad to hear it works. Probably you were simply off to an unfortunate start, so don't worry too much. Good luck with your analysis :)

So, as I mentioned before, that worked beautifully for the human datasets. I'm trying with some mouse data and I'm getting a different error -- I'm not sure if I should make a new question, but I'll post it here: when I run ensembl = useMart("ensembl", dataset = "mmusculus_gene_ensembl") but then when I run mart, I get Object of class Mart: Using the ENSEMBL_MART_ENSEMBL BioMart database, Using the hsapiens_gene_ensembl dataset, instead mmusculus. I tried referencing the archive too, but that didn't work either. Any suggestions appreciated!

turns out one of my input files was null! all good now thank you again.

Those aren't external_gene_name's, they're ensembl_gene_id_versions:
theBM <- getBM(attributes='ensembl_gene_id','hgnc_symbol',
               filters = 'ensembl_gene_id_version',
               values = charg2,
               mart = ensembl)

Note that you'll get more hits if you strip the gene ID versions off:
charg2 = sapply(strsplit(charg, '.', fixed=T), function(x) x[1])
theBM = getBM(attributes='ensembl_gene_id','hgnc_symbol',
              filters = 'ensembl_gene_id',
              values = charg2,
              mart = ensembl)

Okay, thank you! So I tried this -- it's returning theBM (finally!) with the gene names, but when I merge it, the merge still doesn't work :(. Thank you immensely though :)))).

Look at the first few rows of each and it should be apparent why there’s a mismatch.

Hm, maybe merge isn't the right tool. I'm trying to have the HGNC symbols replace my column of Ensembl Gene IDs -- the getBM worked beautifully, but what next?

left_join from dplyr

Okay, please help, how would I implement that?

Sorry for the double comment -- this has been gnawing away at me. I was looking through the biomaRt documentation and I found this line: # rows are sorted by filters, and in different order than the values argument, which makes sense -- my getBM is working. My question is, how do I get the order of getBM to match my results?

Use factors to apply the sorting order from one dataset to another.


text:
What does instantaneous rate matrix mean?
PhyML is a tool used for maximum likelihood estimation analysis.
I was running a PhyML analysis to check the rate of nucleotide change in DNA sequences and in the PhyML output, there are rate parameters and an instantaneous rate matrix. Could you please clarify the difference between the two?

What I've understood so far is that the rate parameters indicate, for ex: that the A -C transition rate is 1.95 and so on. If so what does the instantaneous rate matrix tell us?

Hi @AntonyGuterres did you assign a root/outgroup to this tree before building it?

Reversible At its simplist the GTR rate, is the General Time Reversible model and most importantly infers the matrix is "reversible", thus as many A mutations will go to C as C mutations go to A and the number of mutations is relative to the G-T mutation rate. The G-T mutation rate is one of the least common mutations (transversion) so all other values are above 1. If you plotted this as a matrix the top right-hand matrix is identical to the bottom left hand matrix because its reversible, i.e. A <-> C . This is cool because you don't need to tell the algorithm which sequence is the root.
Directional The "instantaneous" matrix is a directional matrix, thus A -> C is not the same as C -> A. In order to obtain an accurate directional matrix you need to assign a root, its like lets start from the deepest node in the tree (the one between the root/outgroup and all ingroup sequences) and now - if that is correct - we can remove the restrictions of reversibility and impose directionality on the mutation process. This is needed if you are performing molecular clock type calculations. However, the key is assigning the root/outgroup - if it assigns it e.g. by mid-point rooting that is not necessarily a good idea.
Advice. My advice is if you are building a basic tree and haven't assigned a root/outgroup, obviously using GTR which is a good move, then you don't to worry about the intananeous matrix because it might not be meaningful. If you are running a molecular clock/dating/rate estimate type model and a very clear what the root/outgroup is and that is specified in the input, then the instananeous matrix becomes important.

The mutation rates of GTR and instananeous are not directly comparable because the GTR rate is with respect to G <-> T mutation, whilst instananeous is primarily concerned with the directionality of mutation and is a 12-parameter model, as opposed to GTR which 6-perameters, where one of those parameters is the reference.
Thus G-T is reversible in GTR (the reference) and so e.g. A<->C is reversible with respect to G<->T.
In instananeous matrix G -> T and T -> G are two different rates as are A-> C and C -> A, so one approximation based on an assumption with respect to a second approximation based on the same assumption can't directly equate to a mutation rate which removes that assumption.

I'm running a rate estimate model. I have supplied the guide tree for it to follow.
According to the Instantaneous matrix, A->C is 4.38 (and C->A is 0.153) but how is this related to the A<->C of 1.95 which is a GTR rate parameter?

They are not directly comparable answered above because A<->C (GTR) is with respect to G<->T mutation rate (which is assigned to a nominal value of 1).

Okay and how do you obtain A<->C (GTR) from A->C and C->A rates? Is it like a weighted average (weighted wrt base frequency)?

I think it is a separate question, I would do it via molecular clock trees.

Last two clarifications

In the matrix, N->C is high for all N. Does this mean that the a mutation most commonly results in 'C' vs any other base?
Also, a value at (x,y) in the instantaneous matrix represent the probability of base x mutating to base y right?

First question, yes, for the root you have set for the tree (it can change if you used a different root). Second question, yes and likely represent (row, column) - which is important in a directional matrix.


text:
Can't align 2 sequences, (MAFFT killed)
I have 2 .fa sequenes like:
>some description of first...
ACTGACTACGTACGACTACGATCTGACTACACGTAGATAGACTAGTCACTACG
ACTGACTACGTACGACTGACTACATGACTGAC...
> some description of second...
ATCGATAGCTACGATCGATAGTAGATAG....

This file is 461mb (1 chromosome of homo sapiens & 1 chromosome of pan troglodytes) and I want to align it. I am using mafft.
mafft --globalpair --maxiterate 1000 Homo_sapiens_1_Pan_troglodytes_1.fa > ./mafft_results/Homo_sapiens_1_Pan_troglodytes_1_output.fa

And it returns me
generating a scoring matrix for nucleotide (dist=200) ... done
All-to-all alignment.
/usr/bin/mafft: line 2719:  6030 Killed                  "$prefix/tbfast" _ -u $unalignlevel $localparam -C $numthreads $seqtype $model -g $pgexp -f $pggop -Q $spfactor -h $pgaof -A $usenaivepairscore $focusarg _ -+ $iterate -W $minimumweight -V "-"$gopdist -s $unalignlevel $legacygapopt $mergearg $outnum $addarg $add2ndhalfarg -C $numthreadstb $rnaopt $weightopt $treeinopt $treeoutopt $distoutopt $seqtype $model -f "-"$gop -Q $spfactor -h $aof $param_fft $localparam $algopt $treealg $scoreoutarg $focusarg < infile > /dev/null 2>> "$progressfile"

Why? How could I align it for further comparing?
EDIT 1
MAFFT working cool, when I tried to use it on test sequence file
>test1
ATCGATCGATGACTACGATCGACTACGTACG
ACTAGCTAGCTACGCATAGTAGATCGA
ATCGACTAGCTCAGATCGACTAGATCGATCGT
>test2
TTATCCACTGATCGACTCAGATCGACTCAGAC
ATCGATCAGCTACGATGACTACGTAC
ACTGACTAGCTAGACTGACTAGCTACGTCA

and after MAFFT mafft test.fa > output.fa it returned

test1
------atcgatcgat-------------gactacgatcgactacgtacgactagctagc
tacgcatagtagatcgaatcgactagctcagatcgactagatcgatcgt
test2
ttatccactgatcgactcagatcgactcagacatcgatcagctac---------gatgac
tacgtac----------actgactagct-agactgactagctacgtca-

So it works.
I can use all 4 at the example. Now about my large file:
alex@alex:~/Desktop/test/BOTH/homo_sapiens_pan_troglodytes$ mafft --thread 4 --maxiterate 2 Homo_sapiens_1_Pan_troglodytes_1.fa > ./mafft_results/Homo_sapiens_1_Pan_troglodytes_1_output.fa
nthread = 4
nthreadpair = 4
nthreadtb = 4
ppenalty_ex = 0
reallocating...
done.
reallocating...
done.
reallocating...
    done.
reallocating...
done.
reallocating...
...
stacksize: 8192 kb
reallocating...
done.
reallocating...
done.
reallocating...
done.
...
generating a scoring matrix for nucleotide (dist=200) ... done
Gap Penalty = -1.53, +0.00, +0.00

Making a distance matrix ..
    1 / 2 (thread    1)
done.

Constructing a UPGMA tree (efffree=1) ...
    0 / 2
done.

Progressive alignment 1/1...

len1=230480994, len2=223796828, Switching to the memsave mode
STEP     1 / 1 (thread    0)m/usr/bin/mafft: line 2719:  3192 Killed                  "$prefix/disttbfast" -q $npickup -E $cycledisttbfast -V "-"$gopdist -s $unalignlevel $legacygapopt $mergearg -W $tuplesize $termgapopt $outnum $addarg $add2ndhalfarg -C $numthreads-$numthreadstb $memopt $weightopt $treeinopt $treeoutopt $distoutopt $seqtype $model -g $gexp -f "-"$gop -Q $spfactor -h $aof $param_fft $algopt $treealg $scoreoutarg $anchoropt < infile > pre 2>> "$progressfile"

Doesnt matter I do not enter --thread 4 or not, --maxiterate 2 or --maxiterate 1000, it still shows me same error. It takes around 2 minute. Ram usage when running at the first minute is 6656mb, after 1 min 11028mb and then it out of ram and error..
PS
discussion moved here

What machine are you running this on, particularly RAM capacity? Just on a remote off chance are there any "core" dumps, e.g. weirdly large files that suddenly 'appear', especially if they have the name 'core' :-) ?

I am running this on  Ubuntu 20.04, CPU(s): 4, RAM 12 GB, Intel(R) Core(TM) i3-3220.. Do not think that this file could be so large for 12gb ram..

I'd do a full diagnostic and post back, because the first stage is obviously working. E.g. have you successfully got MAFFT working on any other data set?  Are all processors being used? Whats the RAM usage when running? Do any other long processes/calculations suddenly quit (e.g. hardware bug)? How long before the error pops up.  What happens when you greatly lower the maxiterations?

@Michael please see my edit 1. But I still want to compare them

--thread 4 is essential, --maxiterate just try and get it working, not getting it to work optimally. Just make sure there are no other processes running on the machine when --thread 4 is operating, e.g. you are using another machine as your desktop
The gist is as you spread the calculation across your processors (4) and simplify the calculation (lowering the iterations), you are getting further into the calculation, for example now you are getting into gap penalties.
If you are agreed then I would try performing the alignment using just one chromosome of the human genome against the chimp genome and I suspect that would work. Then take it from there.
The objective is to get the calculation working, for example I wouldn't trust a result with max. iterations of 2, but once you have it working you can then start refining the parameters towards optimal results. The i3 processor isn't perfect in context.
Using --thread 4 is super important and you want to increase maxiterate towards the default if the machine supports it.
I suspect the solution is to perform the calculation on an 8+-core machine and set --thread 8 (or higher depending on the number of cores on the machine), at least it will get further into the calculation. With 4-cores and an i3 (diagnostics of RAM usage would be nice) in any case it just doesn't seem to like it. The alternative is to perform this across the genome via sampling the chromosomes on your current i3.

Just to explain why this works. When you have a RAM intensive job and you spread it across cores the RAM overhead massively reduces for each core that is added. Don't ask me how this works in the hardware mecahanics of RAM distribution per core. The reason I know this is because of emperically parallelising Bayesian calculation using each chain on a different core.

I am actually comparing all chromosomes of human vs all chromosomes of chimp. I doesn't clearly understand what do you mean. You chould try it yourself I can wait

Hi @Alexey you have to get it working first. Only when it is working can you move towards the calculation you want. This is not the finished product just a means of verifying it will work. One issue is the i3 processor in the machine. I wouldn't trust a result with max iterations at 2, but once it works you would increase that parameter.

Basically to cut to the chase, I think your machine is underpowered for the calculation you are performing.