text:
Generating multiple NCBI URLs with eutils
I would like to create a URL that takes input from different arrays to generate a search. The data I would like to retrieve is from NCBI. Below you can find the code I wrote. But I can't get it to work; I am stuck. Thank you in advance.
I should have something like this as result but i got nothing. https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=gene&retmode=xml&retmax=10&sort=relevance&term=A20%20AND%20Homo%20sapiens%20%5borgn%5d%20AND%20alive%5bprop%5d
$eutils = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils";
$db = "gene";
$gene_symbols = array("A20","HBB","ABO");
$species = array("homo sapiens", "Escherichia coli","Saccharomyces cerevisiae");
$links = array();
foreach ($gene_symbols as $key => $GS){
    $SP = $species[$key];
    $query = "$GS AND $SP [orgn] AND alive";
    $esearch = "$eutils/esearch.fcgi?db=$db&retmode=xml&retmax=10&sort=relevance&term=";
    #echo $esearch;
    $links[] = get($esearch.$query);
}
foreach($links as $lk){
    xml = simplexml_load_file($lk) or die("feed not loading");
}
var_dump(xml);

As output I got nothing.
But I would like to have this: https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=gene&retmode=xml&retmax=10&sort=relevance&term=A20%20AND%20Homo%20sapiens%20%5borgn%5d%20AND%20alive%5bprop%5d
And this for every gene. For some reason it put them all after each other and the foreach loop doesn't seem to work. I got no errors.
I run this code in php using visual studio code. The result should be visualised in the webbrowser (http://localhost/test.php).

Are you getting no output at all or simply the wrong output?

Please [edit] your question and i) explain how you run this code and on what system; ii) tell us what language it is written in. Is that PHP? iii) Most importantly, explain exactly how it fails. Any errors? Any data at all?

Can you print the links to see if they are valid? That would narrow down where the problem is

As output I get this:http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=bccmtest&retmode=xml&retmax=10&sort=relevance&term=
But I would like to have this: https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=gene&retmode=xml&retmax=10&sort=relevance&term=A20%20AND%20Homo%20sapiens%20%5borgn%5d%20AND%20alive%5bprop%5d
And this for every gene. For some reason it doesn't add the term and the foreach loop doesn't seem to work. I got no errors.
I run this code in php using visual studio code. The result should be visualised in the webbrowser (http://localhost/test.php).

The output you linked to is an error: Invalid db name specified: bccmtest. Your code shows you are setting the db to gene, so what do you get when you run the actual code from your question?

I forgot to change the db name in the link I posted. I would like retrieve the search in the gene database from ncbi. The name of the database is gene.  If I do that i got "Empty term and query_key - nothing to do".

cross posted: https://www.biostars.org/p/311903/

I don't know PHP, so I'm sure there are better ways of doing this, but you had various issues in your code. This seems to be a working version of what you're after:
<?php
$eutils = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils"; $db = "gene"; $gene_symbols = array("A20","HBB","ABO");

$species = array("homo%20sapiens", "Escherichia%20coli","Saccharomyces%20cerevisiae");
$links = array();

foreach ($gene_symbols as $key => $GS){
  $SP = $species[$key];
  $query = "$GS%20AND%20$SP%20[orgn]%20AND%20alive";

  $esearch = "$eutils/esearch.fcgi?db=$db&retmode=xml&retmax=10&sort=relevance&term=";

  echo "$esearch$query\n";
  $links[] = file_get_contents($esearch.$query);
  foreach($links as $lk){
    echo "|------ $lk ----|";
    $xml = simplexml_load_string($lk) or die("feed not loading");
    var_dump($xml);
  }
}
?>

There doesn't seem to be a get() function in PHP, you want file_get_contents. Also, your loops were nested incorrectly and you needed var_dump($xml) to dump the variable $xml not var_dump(xml); which tries to dump the string xml. You're not getting a file here but a string, so you need simplexml_load_string. Finally, you can't have unescaped whitespace in URLs.
Next time, the way to debug this sort of thing is to have it print the URL it is generating and try to visit it manually. You would have seen the spaces and noticed that the working example you show uses %20 instead of space.

I don't know anything about php, but when I copy your code into here it isn't happy about the variable name $query.
This code appears to produce valid urls, but I can't make the get part work:
$eutils = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils"; $db = "gene";$gene_symbols = array("A20","HBB","ABO");
$species = array("homo sapiens", "Escherichia coli","Saccharomyces cerevisiae"); $links = array();

foreach ($gene_symbols as $key => $GS){
   $SP = $species[$key];
   $SearchTerm = "$GS AND $SP [orgn] AND alive";

   $esearch = "$eutils/esearch.fcgi?db=$db&retmode=xml&retmax=10&sort=relevance&term=";

   echo $esearch.$SearchTerm;
}

http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=gene&retmode=xml&retmax=10&sort=relevance&term=A20 AND homo sapiens [orgn] AND alivehttp://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=gene&retmode=xml&retmax=10&sort=relevance&term=HBB AND Escherichia coli [orgn] AND alivehttp://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=gene&retmode=xml&retmax=10&sort=relevance&term=ABO AND Saccharomyces cerevisiae [orgn] AND alive

$eutils = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi";

   $databank = "gene";
   $gene_symbols = array("A20","HBB","ABO");

   $species = array("homo sapiens", "Escherichia coli","Saccharomyces cerevisiae");

foreach($gene_symbols as $value => $GS){
    $Species = $species[$value];

    $esearch = "$GS%AND%alive[prop]%AND+$Species[orgn]";
    $esearch_test = "$eutils?db=$databank&retmode=xml&retmax=10&sort=relevance&term=$esearch";
    $result = file_get_contents($esearch_test);
    echo "Query: " . $esearch_test. "<br \>\n";
    echo "Result: " . $result . "<br \>\n";

    $xml = simplexml_load_string($result);
    echo "Count:<br \>\n";
    echo $xml->Count;
    $C =  $xml->Count;
    echo $C;
}


text:
What are negative FSC-A values in flow cytometry data?
I did flow cytometry analysis of E. coli cells 7 hours post inoculation (1% overnight culture in 3ml).
I am analysing the FCS3.0 files using the flowCore package in R. I observe that FSC-A has negative values in some events, although the corresponding FSC-H values were all positive. These points comprised just a minor fraction of the data (~250/10000 in this case). The FSC voltage was set at 788V and the threshold flow rate was 5000-10000 events/sec.
What do these negative values mean?

FSC-A vs FSC-H: plotted only negative FSC-A values

SSC-H vs FSC-A: red points indicate negative FSC-A values

Why did you delete the SSC-H plot? Also, yould you post the corresponding SSC-A/-W FSC-W plots?

@bukwyrm well, what I claimed in the image that I posted wasn't really true. I don't have the -W data. I tried to see if there is any correlation but the only thing that I could find is that these negative FSC-A points have low SSC-H  but removing these points from the total data does not change the distribution of the latter.

This is probably a matter of calibration. If this is not the raw data (actual voltage curves), look at that.
As an area below 0 is not possible, it has to be a calibration issue, either in the machine, or in the software.
the -H plots are about the maximum height of the voltage spike in the side or forward light detector, while the -W plots are about the width (using some definition of flank-detection, so variation prone).
The area is also computed in various ways: using something like (H*W)/c (with c something like 2) or with an actual computation of the area of the spike, but again this has all kinds of predefinitions in it (what's the beginning/end of a  spike, what's the background-level, ...) Some of these predefinitions you can alter in software, some are hardcoded or even part of the functionality of the hardware, so it might be easy to impossible to get rid of the negative values.
The FSC-H values corresponding to the FSC-A values are in your plot, and positive, so the most obvious solution (a bad algorithm that calculated Area as Height times Width times Magic and did not account for negative Heights (which depending on spike definition, you might get)) is out.

So basically there is some unknown bad algorithm that is responsible for this. Should I consider FSC-H as the reliable estimate and retain these points for further analysis or discard them? As I mentioned in the other comment, removing these points does not change the SSC or FITC distribution.

source: Stack Exchange

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

sort: votes


text:
What are phantom peaks in ChIP-seq?
I've seen two uses of this term that seem to refer to completely different phenomena. Is my understanding correct?

In Jain et al. "Active promoters give rise to false positive ‘Phantom Peaks’ in ChIP-seq experiments" (PMC4538825), a phantom peak is a region on the genome where lots of reads align in a control experiment where there is none of the protein being IP'd, but not in a control experiment without any antibody. The authors (speculatively) attribute this to general stickiness of transcription factors and the transcriptional machinery, which can bind with potentially any antibody and get IP'd.
The ENCODE quality metrics (link) say "a cross-correlation phantom-peak is ... observed at a strand-shift equal to the read length." They are referring to a high correlation of x[i] and y[i + r] where r is the read length, x[i] contains counts of tags aligning to the Watson strand at coordinate i, and y[i+r] contains counts of tags aligning to the Crick strand at coordinate i+r. They attribute this phenomenon to variations in mappability. (I don't understand their explanation at all, but that's a separate question.)

This is almost right, but it's missing a little context. In ChIP-seq, it's typical to shift the positive and negative reads because otherwise the peaks end up offset from each other. The amount of shift is determined by trying different shifts and measuring how well the read densities correlate between positive and negative strands. In the question, the first sense talks about a peak on the plot of read density vs genomic position. The second has different x and y axes because it's talking about the cross-correlation between strands. X is the offset between positive and negative strands, and y is the correlation observed at that offset. Look at figures 4D,E in this ENCODE paper for visuals.
Landt, S. G., Marinov, G. K., Kundaje, A., Kheradpour, P., Pauli, F., Batzoglou, S., ... & Chen, Y. (2012). ChIP-seq guidelines and practices of the ENCODE and modENCODE consortia. Genome research, 22(9), 1813-1831.
https://genome.cshlp.org/content/22/9/1813.full.pdf

i'm not sure what the best terms are for strands -- "positive"/"negative"; "watson"/"crick"?

I think positive / negative (as you have here) is reasonable, especially as it relates to gene-relative direction.

Positive/negative, sense/antisense, top/bottom, forward/reverse, plus/minus... All describe the same thing.

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

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


text:
Merge 2 VCFs from different variant callers
I am working with WES data for detection of somatic variants and I have used two variant callers because no variant caller is complete in itself. I have used GATK Haplotypecaller for small variants like snp,ins,del etc. And for large variants I have used pindel (specially for SVs detection). Now I have two vcf file for each patients, one from gatk and another from pindel. Although I know various vcf-merging tools but is it right to merge two different vcf file from different callers. Do they give right and proper merged vcf file. Are those vcf-merging tools SV-aware and consider only trivial call for the same variant represented differently in different callers. I know these are too many questions for a single thread but all these are related to each other. I am in dilemma whether I should merge them or use them separately.
Thanks.

The merge utility from SURVIVOR seems to be what you are looking for. I did not try it yet, but it seems to be especially designed to handle SV.
The simplest thing to do would probably to first merge your patients files together for each tool and have one multi-patient file per variant caller. Then I would merge the multi-sample files of the different tools.
Here is a detailed blog post describing the different challenges associated with merging structural variants from different callers and explaining the steps they took to address it. All the code to reproduce the pipeline is available on a github repo.

Thanks for your reply and sorry for late response. Survivor is for merging SVs from different variant caller but can we merge small and large variants vcf file using this?

I am not sure, but It would seem it handles both small and large variants by looking at the source code. They also have a wiki detailing the methods. I think you might want to try on a test file and see if it does what you expect.

Thanks. I'll try and let you know if it works well with both variants or not.

Hey!! I know it is very late to respond but I was quite busy in other stuffs. I tried this tool and found that this tool extract SVs from both VCF file and merge them in a single file. But I need merging tool which combine all variants from both VCF file including SVs and merge than into single VCF file not just SVs. Unfortunately, that does not work. Any other idea... Thanks


text:
Get gene names corresponding to CPG probes
I have downloaded KIRC Methylation data (450K) using TCGA toolbox. The beta values given are with reference to CPG probes. Is there any package or function in R through which I can get the gene names corresponding to these probes? I was able to get the gene I.Ds using getMappedEntrezIDs. Is there any function to map these I.Ds to respective gene names?

Here's a Python script (e.g., called map_entrez_to_ensembl.py that uses mygene.info to convert a list of Entrez IDs to Ensembl gene names:
#!/usr/bin/env python

import sys
import mygene

mg = mygene.MyGeneInfo()

genes = ['3815',
         '3816',
         '2341']

results = mg.querymany(genes, scopes=['entrezgene'], fields=['ensembl'], species='human', verbose=False)

for res in results:
    entrez_id = res['query']
    if 'ensembl' in res and 'gene' in res['ensembl']:
        sys.stdout.write('%s\t%s\n' % (entrez_id, res['ensembl']['gene']))

Result:
$ ./map_entrez_to_ensembl.py
3815    ENSG00000157404
3816    ENSG00000167748

Looks like there's an R-based interface to mygene.info, if you really need to use R.
You could use their API test page to see how data are structured and what keywords are available for different kinds of searches.
For instance, if you wanted aliases, in addition to Ensembl names:
#!/usr/bin/env python

import sys
import mygene

mg = mygene.MyGeneInfo()

genes = ['3815',
         '3816',
         '2341']

results = mg.querymany(genes, scopes=['entrezgene'], fields=['ensembl', 'alias'], species='human', verbose=False)

for res in results:
    entrez_id = res['query']
    if 'ensembl' in res and 'gene' in res['ensembl'] and 'alias' in res:
        for alias in res['alias']:
            sys.stdout.write('%s\t%s\t%s\n' % (entrez_id, alias, res['ensembl']['gene']))

Then:
$ ./map_entrez_to_alias_and_ensembl.py
3815    C-Kit   ENSG00000157404
3815    CD117   ENSG00000157404
3815    PBT     ENSG00000157404
3815    SCFR    ENSG00000157404
3816    KLKR    ENSG00000167748
3816    Klk6    ENSG00000167748
3816    hK1     ENSG00000167748


text:
How to add data to AnnotationHub
I recently started using the Bioconductor AnnotationHub package for accessing lots of common datasets, but I've noticed that there are a lot of newer databases that have not made it into AnnotationHub's collection.
Has anyone had any luck getting databases added to AnnotationHub or adding them manually yourself (I guess via the package process?)?  Is it painful or relatively straightforward?
I love the annotationhub interface, but I'm wondering if it is way too inefficient to add new resources rather than just manually downloading and curating myself.  Thanks!

As the description of the package says:

The AnnotationHub web resource provides a central location where
  genomic files (e.g., VCF, bed, wig) and other resources from standard
  locations (e.g., UCSC, Ensembl) can be discovered.

If the new resources aren't in the standard locations it won't be listed by AnnotationHub. The easy/common way to have the data accessible through AnnotationHub is make it available through these standard locations (before publications if you want to give early access).
Another method would be to follow this steps of the vignette to add new resources, but it might involve creating a new package.

Yeah, the weird thing is that the data I'm speaking of is the BLUEPRINT epigenetic datasets from hematopoietic cells, and this data is already in Ensembl, so I thought it would be accessible by AnnotationHub, but it is not yet.  Looks like I'll just have to do it manually until they update it on AnnotationHub. Thanks for your insight!

@Reilstein If the addition is recent maybe it hasn't been made it yet to the Hub (I don't really know how does it work internally). You could also contact the maintainers to know more about this issue, perhaps is there an issue they aren't aware of.


text:
How to predict the distance between two residues in a protein sequence
I have some protein sequences and I need to predict the physical distance between each two residues in the protein when folded correctly in 3 dimensions; I need to predict this distances just by having the protein sequence and not the structure.
I have a group of pdb IDs and I have got homologs for them and made a multiple sequence alignment for each query sequence; But later I need to imagine that I don't have any structure for being able to predict the distance between residues in query sequences without any known structure.
There are some articles for solving this issue using machine learning methods but I don't want to train with a large numbers of structures. Any one knows any way to solve this problem by using mathematical or statistical methods?

Are there any similar sequences with known structure?

Please add that detail into your question. Also, pelase clarify that you are referring to the physical distance between the residues in the protein when folded correctly in 3 dimensions (at first I thought you meant in the linear sequence).

Between ALL pairs of residue-residue combinations? You can't use statistics without a large enough dataset and expect reliable results... BTW, I might be outdated but I don't remember having seen de novo folding success after Brian Kuhlman's, 2003 science paper with David Baker. If you only need to compute the distance within short fragments, you could estimate the secondary structure with psipred and use this information to estimate, for example, an average distance between residue i and residue i+4

@aerijman I think you are mixing up de novo protein design and de novo protein structure prediction. Design of novel folds hasn't had much success as you say, but structure prediction of novel folds is getting to a decent level now. See the CASP website for some nice examples.

Developments in this area over the last year are worth pointing out. All use residue-residue covariation information from a protein family so require a multiple sequence alignment to work, though it sounds like you have that.
The most accurate way to do this would be to use pre-trained machine learning methods, which would not require you to train your own version. Examples include:

DMPfold from our lab, which is freely-available open source.
The RaptorX method from the Xu lab.
The state of the art is AlphaFold from DeepMind, though this has not been published yet.

These predict the probabilities of a residue-residue distance being in various distance 'bins', e.g. the probability of being 5 - 6 Angstrom apart. They don't predict beyond ~20 Angstrom.
If you specifically wanted a statistical, non-trained method, then consider:

PSICOV from our lab.
CCMPred.
EVfold.

All of which are available in some form. The problem here is that they all predict binary contacts at 8 Angstrom, whereas it sounds like you need more fine-grained distances.
In a general sense, you can use any 3D structure prediction method and then read the distances back from the predicted structure.


text:
10X cellranger count, [error] The chemistry was unable to be automatically determined
While running cellranger count using the following .slurm file:
#SBATCH -p standard
#SBATCH -A overall
#SBATCH --time=48:00:00
#SBATCH --mail-type=end
#SBATCH --mail-user=nv4e@virginia.edu
#SBATCH --output=generate_count_%j.out
#SBATCH --error=generate_count_%j.err
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=20
#SBATCH --mem=120000

# Run driver program
FASTQS="/scratch/nv4e/kipnis/AHN3LLBGX5/outs/fastq_path/HN3LLBGX5"
TR="/scratch/nv4e/kipnis/refdata-cellranger-mm10-1.2.0"

cellranger count --id=17R_1 --transcriptome=${TR} --fastqs=${FASTQS} --sample=17R --expect-cells=2000

I got the error:
>2018-04-28 20:04:37 [runtime] (progress)
ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMI
STRY.fork0.chnk0: Indexing genome...
2018-04-28 20:04:49 [runtime] (progress)        ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMISTRY.fork0.chnk0: Building transcriptome...
2018-04-28 20:04:58 [runtime] (progress)        ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMISTRY.fork0.chnk0: Building kmer index...
2018-04-28 20:05:25 [runtime] (progress)        ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMISTRY.fork0.chnk0: Mapping reads...
2018-04-28 20:05:28 [runtime] (chunks_complete) ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMISTRY
2018-04-28 20:05:28 [runtime] (run:local)       ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMISTRY.fork0.join
2018-04-28 20:05:31 [runtime] (join_complete)   ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMISTRY

2018-04-28 20:05:34 [runtime] (ready)           ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.SETUP_CHUNKS
2018-04-28 20:05:37 [runtime] (split_complete)  ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.SETUP_CHUNKS
2018-04-28 20:05:37 [runtime] (run:local)       ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.SETUP_CHUNKS.fork0.chnk0.main
2018-04-28 20:05:40 [runtime] (failed)
ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.SETUP_CHUNKS

>[error] The chemistry was unable to be automatically determined. This can happen if not enough reads originate from the given reference. Please verify your choice of reference or explicitly specify the chemistry via the --chemistry argument

What is going wrong here? How could I correct it?

You presumably aligned against the wrong reference transcriptome.

It is very weird: the thing that the first time I ran it, it worked. Then, because we got bad results after running cellranger count ('bad' means biologically not what we want to see. It would mean that the biology was done in the wrong way), we transferred the raw data again, generated fastq files again and somehow now I am having this issue for the first replica, but not for the second one which ran successfully. I am wondering where in this chain of actions I did a mistake. Can it be that the data got corrupted while I was transferring it, however, fastq generation worked still?

The issue was that I supplied wrong indexes during fastq files generation, however, they were generated successfully, and only the next step failed with not that informative error.


text:
What does the setSize ordering mean?
I am running gseGO function from clusterProfiler and getting an output where terms are ordered by the setSize parameter. What does this ranking mean? Are terms ordered by the sizes of the gene sets? If yes, then what is the ultimate value of such ranking? My goal is to get the most relevant terms, but I am not completely sure that such an ordering based on the setSize is actually the right one.

Did you read the clusterProfiler vignette and messages of the author in the support site and its own blog?

Sure, neither of them provide enough info.

If you post the head of your output you are talking about, it could be easier to clarify the doubts

As you can see on this blog post of the author, the result is ordered by pvalue, but it can coincidentally sort by setSize. The setSize means that that term has so many genes on under that category.

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


text:
Finding new drug-protein targets/interactions using a known one as base
Given a known interaction of a molecule with some proteases (having X-rays images too) is it possible to predict interactions with other proteins either using their 3D structure or not?
If anyone has already tried such techniques I would like to hear a little brief on that subject. Which are the possible techniques, computational power needed and so on.

If you have information on which part of the protein are interacting, and its sequence, you can search for homologous or similar sequences in other proteins. This would give you a list of potential targets and probably a similarity score which you can also use. Then you would have something to run simulations on.
And I think that with 3D structure you should be able to run simulations in silico in this case (eventhough it is probably more common to look for the ligands to the protein, not proteins to the ligand).
In silico simulations are computationally intensive, but I´ve also seen some being run on a laptop, so you don´t necessarily need access to HPC-clusters, or supercomputer.

This sounds like something you could possibly do, but is it safe to assume that the evolution of interactions is way faster than the sequence evolution? Is there any described method or paper using the approach you suggested?

@Kamil S Jaron And why we want to assume that?

You can find information in uniprot.org database (on structure, interactions, domains and motifs and similarity) and it also have supporting data on the data sources and literature citations.


text:
Will blastn remove sequences from a search with low identity?
I'm using this command (excuse the duplicate naming, I know it's bad form):
blastn -query cm-seqs/combined_seqs.fna -db combined_seqs.fna -out cm-matched.txt -num_alignments 1 -outfmt 10

to take a set of fasta sequences to blast against another data base I built with blast.
My query file has 4,364,417 sequences. The files resulting from the above command spat out 4,362,639 results. The difference is about 1800. Is this due to poor alignments for some sequences that didn't pass any of blasts default parameters? I assumed the result would be 1 "best" match in the database for each query sequence.

-num_alignments 1 is not appropriate for that output format i think, try -max_target_seqs 1

Yes, if a query sequence has no good alignment, blast will not return a result for it. Your command will indeed output only the best hit for each query sequence, but there is no reason to assume all of your sequences will have a match.
Those that don't have a match won't be included in the output. This sounds perfectly normal unless you are blasting the same input file against a DB created from that file. Otherwise, you don't expect all of your sequences to necessarily have a match.
Also note that since the default blastn settings are quite permissive (usually an e-value of 10, for example), don't make the mistake of assuming that a blast hit means anything. The details will always depend on the specific search you are running (I have discussed this a little in my answer here), but in the vast majority of cases, any hit with an e-value > 1 can safely be ignored. In most cases, you are only interested in e-values well under 1, actually. So you probably already have all sorts of spurious results in there anyway.
For more details, please edit your question and explain what the query and DB files are.

Great! That was a perfect response.

The databases I created were the bacterial 16S metagenome of an additive, and the 16S metagenomic community of our experimental subject during the experiment. The microbiome of the subjects were blasted against the additive to track what in the additive was able to persist in a foreign environment.


text:
Human SNP databases other than dbSNP
I'm working on a project to identify variants as being either germline or somatic based on nearby SNPs, and I was wondering if there were any databases of SNPs for Homo sapiens other than dbSNP? I'm not looking for databases that contain additional information (pathogenicity etc.) on SNPs that are already in dbSNP, instead I want to know if there are any other databases that contain SNPs that aren't found in dbSNP.
I've been searching online and so far the only databases I can find are either deprecated or integrated into dbSNP, so any suggestions would be great.

I'm pretty sure there are no public databases of variants that aren't included in dbSNP. That is the variant database so that's the one you want to use.
While there are various other databases around, they tend to use dbSNP rsIDs as well so any variants in those will also be in dbSNP. So while dbSNP is nowhere near complete, of course, it's still the most complete one you're likely to find.

Fair enough, just thought I'd check. Thanks!

Just adding to that, some of the more specialised databases (including OpenSNP) can be found for example here: https://biokeanos.com/search?q=SNPs+snp+%2Bhuman.


text:
Declaring an array of strings and integers
I have two csv file - one is a list of strings and other has both strings and floating values. If the strings in both the files match I need to copy the complete row from the second file to a temporary array. I have declared a matrix with n rows and m columns and saved is as dataframe temp <- as.data.frame(matrix(nrow = n, ncol = m)). I am using for loop to iterate over all rows and strcmp to compare the strings. But when the data is written to temp the strings are stored as integers. How should I declare this array? Is it possible to write this loop without declaring an array in R?

This question is unrelated to bioinformatics. You should ask somewhere more appropriate (e.g. stackoverflow). You might also want to add a minimal working example, it will be easier for other people to help you (see: https://stackoverflow.com/help/mcve). Good luck !

I'm downvoting because this isn't really a bioinformatics question, but I suggest you take a look at dplyr::semi_join

The merge function of R might be useful. The function expects two data frames, with [optionally] columns to compare, and it will output all fields from both data frames where the comparison columns match.
> a <- data.frame(test=1:10);
> b <- data.frame(test=c(sample(1:10,5),sample(c("A","B","C","D","E"))),
                  output=rnorm(10));
> merge(a, b, by="test");
  test     output
1    1 -1.3800443
2    2  0.2934274
3    3  0.0868700
4    9  1.1766914
5   10  0.3793448

The merge function can be customised quite a bit (e.g. if you want NA set for non-matching rows), see ?merge for more details.

When you read in your two CSV files, you can explicitly declare the variable types for each column:
> a <- read.csv(..., colClasses=c("character"), ...)
> b <- read.csv(..., colClasses=c("character", "numeric", "numeric", ...), ...)

You could also investigate as.is or stringsAsFactors, although this can cause problems where there are duplicate strings. Using colClasses is more explicit.
Read the help documentation in ?read.csv for a more complete explanation of these operators.
You can filter with the %in% operator in match:
> strings_of_interest <- as.vector(a$V1)
> b.filtered <- b[b$V1 %in% strings_of_interest,]

If I am reading your question correctly, perhaps you can then take b.filtered and turn that into the matrix you want.
Avoid use of for loops where you can. They can slow R down considerably. A good read on this subject is in the third chapter of Patrick Burns' The R Inferno.

Is there a reason for the downvoting?

I think it is bad for a community when downvoting without any explanation. Please explain yourself.


text:
Convert CD (cluster of differentiation) into gene symbol
I want to convert names according to the "cluster of differentiation" nomenclature into accepted gene symbols. For example, CD3 should map to "CD3a;CD3b;CD3c;CD247" and "CD45RO" should map to "PTPRC". If the surface antigen isn't coded for by a gene, it should map to "" or something similar.
I've seen Obtain CD (cluster of differentiation) marker expression information, but the Abcam table used there does not give the gene symbol -- it gives a list of alternate names. I've seen table 1 of this paper, but it contains errors -- for example, TNFRSF11 should be TNFRSF11A. I could copy this table by hand, or I could download this one, though that's annoying because it doesn't have one row per CD.
Is there a canonical (or just a convenient and reliable) way to do this?

I ended up using option 4 and parsing the columns for "Synonym" and "Previous Symbol".
https://www.genenames.org/cgi-bin/genefamilies/set/471


text:
Keep Format and Individual fields when annotating VCF with VEP
I'm currently updating my Variant Calling Pipeline by switching the VCF annotating software from Annovar to VEP for a variety of regions, not least how easy it is to annotate with HGVS notation and keep datasets up to date in VEP.
For the most part everything is running smoothly, with the exception that some of the data in the VCFs is lost during annotation (and conversion to tsv). The VCFs are created with GATK's UnifiedGenotyper and include a 'Format' column where each value is 'GT:AD:DP:GQ:PL' and a column named after the Individual, which contains semicolon-separated data that corresponds to the Format column (i.e. Genotype;Allele Depth;Depth;Genotype Quality;Phred-likelihood). When I annotate with VEP none of this data is carried over to the output file as it would be in Annovar, leaving me with an annotated file that has no information on read depth, genotype or any of the other data in the two lost columns.
I've included the command I'm currently using for annotation:
./vep -i RM0108.vcf --cache --force_overwrite --tab --merged --variant_class --sift b --polyphen b --hgvs --symbol --canonical --check_existing --af_1kg --af_gnomad --humdiv --pick -o RM0108.tsv

I can't find information on this in the VEP documentation or elsewhere online. I could write something to take the relevant information from the VCF and add it to the tsv after VEP has finished running, but it seems like there may be an easier solution that I'm missing, so any help would be appreciated.
I've also posted this question on BioStars Link Here

The VEP TSV format can only keep in its own specified columns. If you want to maintain the data from your original input, get your output in VCF. It will add the VEP annotation to the INFO column, and keep all the stuff you already have there.


text:
Samtools/bcftools calling indels in noisy reads
I have very noisy nanopore reads and am trying to call SNPs/indels.
I'm runinig into some trouble when truing to use the samtools mpileup | bcftools call combination.
It seems to incorrectly be calling very long indels where it looks like there is no support. For example:
samtools mpileup -f ref.fa sample.bam -r Chromosome:198940-198940 produces:
Chromosome    198940    C    37    .,.-1A.-1A...-1A.-1A.-1A.,-1a.-1A.,-1a,-1a.,.-1A,,.-1A,,+1a,+1a,.-1A,,-1a.,-1a,-1a,,-1a,-1a.-1A,-1a.    06?5/<:8>/4<?3691465<3354.08:23:11@3C

This looks fine... however when I try to call the genotype with bcftools I get a very long indel.
samtools mpileup -ugf ref.fa sample.bam -r Chromosome:198940-198940 -t AD | bcftools call -mv produces:
Chromosome    198940    .    CAAAAGACGTCATCGACGTACGGTCGGTGACCACCGAGATCAACACGTTGTTCCAGACGCTCACCTCGATCGCCGA    C    225    .    INDEL;IDV=32;IMF=0.5;DP=62;VDB=0.145385;SGB=-0.662043;MQSB=1;MQ0F=0;AC=2;AN=2;DP4=0,0,7,2;MQ=60    GT:PL:AD    1/1:255,27,0:0,9

Using htsbox pileup -f ref.fa sample.bam | grep 198940 it shows there is one read with this indels together with many other reads with smaller indels.
Chromosome  198940  C   C,C-75AAAAGACGTCATCGACGTACGGTCGGTGACCACCGAGATCAACACGTTGTTCCAGACGCTCACCTCGATCGCCGA,C-4AAAA,C-2AA,C-1A,T,C+1A,C+1T,C+3GGT 0/4:25,1,1,2,24,5,2,1,1

Is there any reason why it would be grouping these reads together into the big indel?
Any idea what is going wrong here?

It's likely that the caller is being overly sensitive to deletions, which are the most common error mode for nanopore reads, even three years later in 2021.
What I'd guess is happening is that the caller is seeing a string of consecutive 1bp deletions at high frequency, and treating these as a consensus deletion covering many more bases. Looking at the pileup calls for the surrounding bases may help to diagnose this.
My general advice to people doing variant calling with nanopore is to be suspicious of deletions, especially if they include homopolymers. That usually means trying to confirm any important variants with another technology.
Newer, more accurate, base calling models are being created all the time by ONT, so if you still have the raw FAST5 data, it might be worth it to run through the most recent base caller.


text:
Finding proteins in a sequence
I have a file with protein sequences, all on a single line, like this:
>Contig_1 1FFPMPQ*CCLFSSSYLPA*AIMVTGSCQTLVQDKV*CLASCHVMIPPCVFHPIYPLFTWSPLHPIISGCF*AVCMCSWIHALFVEIILSWSLSLPFTFEISSD>Contig_2 1FFPCHNSVVYFLLLICLLKPLWSQGVVRHWYRTRCDAWLAVML*YLHVSFIPSIHFLLGRLCTLLYLDVFELYACAAGFMPCLLR*FSLGHSLFLLPLKFPATVKFFSLLGQIMNNISLMNDLQTGEKLQETH*EVEDKSALLVSIIESRCCVTLQPQDFSN*HLGKL*KSVGQKLK*VNDQITY*

I want to print all sequences between an M residue and a * that are at least 6 residues long, and I would like the output like this:
>contig_1.1.1
M.......*
>contig_1.1.2
M........*
>contig_2.1.1
M............*
>contig_2.1.2
M......*

Where the first number after the underscore is the number found in as the first character of the sequence (it is the reading frame) and the second number is the sequence number. and the last number is the fragment number. So, since the 1st sequence starts with;
>Contig_1 1FFPMPQ

The name will be >Contig_1.1 and each fragment found will be >Contig_1.1.1, >Contig_1.1.2, ... , >Contig_1.1.N.
After implementing the suggestion in the answer below and escaping the * in my regular expression, my current code is:
import sys
import re

myprotein = ">Contig_1 1FFPMPQ*CCLFSSSYLPA*AIMVTGSCQTLVQDKV*CLASCHVMIPPCVFHPIYPLFTWSPLHPIISGCF*AVCMCSWIHALFVEIILSWSLSLPFTFEISSD>Contig_2 1FFPCHNSVVYFLLLICLLKPLWSQGVVRHWYRTRCDAWLAVML*YLHVSFIPSIHFLLGRLCTLLYLDVFELYACAAGFMPCLLR*FSLGHSLFLLPLKFPATVKFFSLLGQIMNNISLMNDLQTGEKLQETH*EVEDKSALLVSIIESRCCVTLQPQDFSN*HLGKL*KSVGQKLK*VNDQITY*"
# parse record
my_protein = []
for line in myprotein:
    if not line.startswith('>'):
        my_protein.append(line)
    else:
        header = line

# build sequence
my_pro = ''.join(my_protein)

# cleave sequence
peptides = re.findall(r'(M[A-Z]{6,}\*)', my_pro)

# build a formatted list
records = ["%sNo.%d\n%s\n" % (header, idx + 1, fragment) for idx, fragment in     enumerate(peptides)]

# print list
for record in records:
    sys.stdout.write("%s" % (record))

And it prints:
 >No.1
 MVTGSCQTLVQDKV*
 >No.2
 MIPPCVFHPIYPLFTWSPLHPIISGCF*
 >No.3
 MNNISLMNDLQTGEKLQETH*

How can I change this to get my desired output format?

re.findall('M[A-Z]{6,}*', seq)

My apologies. the back-slash in my code got stripped for some reason. Try re.findall('M[A-Z]{6,}\*', seq)

Hi and welcome to the site! When posting this sort of question, please always include the errors you get. Telling us "this doesn't work" is really not very helpful, it is always better to explain how it fails.

Sorry for posting like that I just update my question so it would be great if you could help. Thanks

The * character has a special meaning in regular expressions, it means "0 or more". So when you write {6,}*, the regex engine tries to apply the "0 or more" to the {6,*} and gets confused:
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python3.6/re.py", line 222, in findall
    return _compile(pattern, flags).findall(string)
  File "/usr/lib/python3.6/re.py", line 301, in _compile
    p = sre_compile.compile(pattern, flags)
  File "/usr/lib/python3.6/sre_compile.py", line 562, in compile
    p = sre_parse.parse(p, flags)
  File "/usr/lib/python3.6/sre_parse.py", line 855, in parse
    p = _parse_sub(source, pattern, flags & SRE_FLAG_VERBOSE, 0)
  File "/usr/lib/python3.6/sre_parse.py", line 416, in _parse_sub
    not nested and not items))
  File "/usr/lib/python3.6/sre_parse.py", line 619, in _parse
    source.tell() - here + len(this))
sre_constants.error: multiple repeat at position 10

In order to use special characters literally, as simple strings, you need to escape them by adding a \ before them. So try:
import re
myprotein = ">Contig_1 1FFPMPQ*CCLFSSSYLPA*AIMVTGSCQTLVQDKV*CLASCHVMIPPCVFHPIYPLFTWSPLHPIISGCF*AVCMCSWIHALFVEIILSWSLSLPFTFEISSD"
protein = re.findall('M[A-Z]{6,}\*',myprotein)
print(protein)
['MVTGSCQTLVQDKV*', 'MIPPCVFHPIYPLFTWSPLHPIISGCF*']

To get the output you asked for while reading from a file, you can do:
import sys
import re

proteinFile = sys.argv[1]

with open(proteinFile, 'r') as f:
    for line in f:
        ## Split the long input line into the list 'sequences'
        sequences = line.rstrip('\n').split('>')

        for i, k in list(enumerate(sequences)):
            ## Skip the first element which will be empty.
            ## there must be a pythonic way of doing this, but I'm
            ## a Perl hacker, so I don't know it.
            if len(k) == 0:
                continue

            seqName, seq = k.split(' ')
            ## Get the frame
            frame = seq[0]
            ## Remove the frame from the sequence
            seq = seq[1:]
            # cleave sequence
            peptides = re.findall(r'(M[A-Z]{6,}\*)', seq)

            for idx, peptide in enumerate(peptides):
                print(">%s.%s.%d\n%s" % (seqName, frame, idx+1, peptide) )


text:
Is there an efficient way to check an input BAM in R?
I'm writing a function in R for an R package which takes as input a BAM.
my_func = function(my_bam){
    if (!file.exists(my_bam)){
        stop("Input a valid BAM file.")
    }

    ## do awesome stuff with BAM

}

Trying to write robust software, I would like to do more of a check whether this BAM is actually a BAM. Above, I'm just using file.exists().
I could also do something simple like "check whether the file has the extension *.bam" with
grepl('*.bam', my_bam)  ## TRUE for file.bam, FALSE for file.txt

However, is there a quick way to check that the input is actually a BAM? I would prefer to throw an error here, if the user inputs a *txt or *bedpe.
I don't think it's necessary to validate the BAM per se, whereby I believe most users rely upon Picard:
https://broadinstitute.github.io/picard/command-line-overview.html#ValidateSamFile
My first thought was Rsamtools for this operation:
https://bioconductor.org/packages/release/bioc/html/Rsamtools.html
However, it doesn't appear that there is a function here which checks the user inputs a BAM, e.g. BamFile() doesn't appear to throw an error
> library(Rsamtools)
> fake_bam = 'nonsense_string.bam'   ## a string I made up, not a real BAM
> BamFile(fake_bam)
class: BamFile
path: nonsense_string.bam
index: NA
isOpen: FALSE
yieldSize: NA
obeyQname: FALSE
asMates: FALSE
qnamePrefixEnd: NA
qnameSuffixStart: NA

How could I "check" that the user is inputting a BAM file? Perhaps pipe to samtools (which doesn't feel very stable in terms of software)?

You won’t know whether it’s a valid BAM file until you attempt reading it in full; even Devon’s method (which I recommend) only does a superficial check — i.e. it checks whether (parts of) the BAM header can be read correctly. However, this doesn’t guarantee that the file isn’t truncated or contains invalid data.
This is of course a more general problem, but it’s perhaps a bit stranger for BAMs (and other biological formats), as it’s relatively frequent to end up with technically invalid files that were produced by existing software. The only recourse is actually trying to use the file, and to handle errors appropriately.
For a very quick check™ you could also pass the file through gunzip and test whether the “magic bytes” at the beginning of the file equal 'BAM\1':
bgzf_stream = gzfile(filename, 'r')
on.exit(close(bgzf_stream))
magic = readChar(bgzf_stream, 4)
stopifnot(identical(magic, 'BAM\1'))

Take your BamFile(fake_bam) and pass it through open(). That will return a value if fake_bam is a real BAM (or CRAM or maybe SAM) file and through an error otherwise. You can then wrap all of that in a trycatch():
trycatch(open(BamFile(fake_bam)), do something, do something else)

Alternatively, ditch the trycatch() and let open() return the error if an invalid file in entered.

Just expanding upon the answers by @DevonRyan and @KonradRudolph, here's one way to put this in an R package:
Either using open(BamFile() from Devon's answer, which needs Rsamtools:
my_func = function(my_bam){

    out = tryCatch(open(Rsamtools::BamFile(my_bam)), error = function(e) e)
    if (inherits(out, "error")){
        stop("Input a valid BAM file.")
    }

    ## do awesome stuff with BAM
}

Or from Konrad's answer
my_func = function(my_bam){

    magic = readChar(gzfile(filename, 'r'), 4)
    if (!identical(magic, 'BAM\1')){
        stop("Input a valid BAM file.")
    }

    ## do awesome stuff with BAM
}

R has a function that specifically performs any(class(out) == 'error'): inherits(out, 'error'). Besides being more readable this is also potentially more efficient (though it won’t matter here).

@KonradRudolph You're right, I wasn't thinking clear. Edited.


text:
FeaturePlot from Seurat: change its title
FeaturePlot is a function in Seurat package.
And in the vignette it is written that if we specify parameter do.return = TRUE it should return ggplot2 object. It is not working. My goal here is just to change the title of the plot. In case of violin plot I can do the following:
VlnPlot(object = seurat_object, features.plot = id, do.return = TRUE) + labs(title = endothelial_symbols[1])

But with FeaturePlot similar code fails to work:
FeaturePlot(object = seurat_object, features.plot = id, cols.use = c("grey", "blue"),
                          reduction.use = "tsne", do.return = TRUE) + labs(title = endothelial_symbols[1])

Giving the error:

Error in FeaturePlot(object = seurat_object, features.plot = id, cols.use = c("grey",  :
    non-numeric argument to binary operator

Any suggestions would be greatly appreciated.

I tried with some data that I have and this is working for me:
p <- FeaturePlot(object = seurat_object, features.plot = id, cols.use = c("grey", "blue"),
                      reduction.use = "tsne", do.return = TRUE)

lapply(p, function(x){x + labs(title = endothelial_symbols[1])})

I think it is because FeaturePlot returns several ggplot objects according to the ids you put. This is returned as a list of plots so you have to iterate each one for labeling.

This can be “simplified” to lapply(p, `+`, labs(title = endothelial_symbols[1])).


text:
Hierarchical models with limma?
I have a dataset with (microarray) gene expression data that was sampled from the same individuals at multiple timepoints. Our exposure is a continuous variable, and because this was an observational study there is no consistent pattern of change in exposure over time (eg. one person's exposure may continuously increase, another's may go up and down, another may have no changes, etc.). The exposure is not expected to have any cumulative effect, so different timepoints within a person do not influence each-other. Our goal is to find genes where expression for the gene is associated with the level of exposure.
We are thinking of analyzing this as a hierarchical model, with timepoints nested within person (so a mixed model with a random effect for person). We will also need to control for a few covariates, some of which will be the same at all timepoints (race, sex) and some of which may change between timepoints (use of particular medications).
Is this something that it is possible to use limma for? I found that limma has a function called duplicateCorrelation() intended for block designs that some people have tried to use to fit mixed models, but I am not sure whether it can be used for this. Or if not limma, then what would you suggest?

Yes, you can use limma for this mixed model approach. Like you suggest, the random effect (persons) can be put in duplicateCorrelation().
Here is a similar example with RNAseq data, on bioconductor support site.

From the help page of duplicateCorrelationEstimate the correlation between duplicate spots (regularly spaced      replicate spots on the same array) or between technical replicates      from a series of arrays.. So I wouldn't advise to use to account for this effect with this function.

@Llopis, did you also read Gordon Smyth's answer in the bioconductor link? Gordon Smyth is the creator of limma, and is professor in gene expression statistics, he knows what he is talking about.

@Llopis It's not explicitly stated as such, but it's another use for that function (see, for example https://support.bioconductor.org/p/59700/ ).

@Devon, I know but it makes it harder to understand what is going on on my POV it would be easier to have it in the design matrix. b.nota, I did and I found that he also says that "You can't test for factor A. It doesn't really make sense to test for a random effect." Where A was the random effect on that question.

@Llopis, putting everything in a design matrix won't make it a mixed model but a linear model with only fixed effects. In the question in the link, the question was to get the effect of A, B, and C, and Gordon explains that it doesn't make sense to get the effect of A (random effect).

Ok, then I need to review my understanding of mixed models, because I thought that adding them to the design matrix was a good way to handle them :\

Don't worry @Llopis, most of the times you will have to deal with fixed effects. In OP's example here, they observed that the effect is random (up, down, or no change), that's when you have to consider using mixed models.


text:
Off-by-one error when mapping Ensembl IDs to gene symbols
I have a matrix with ensemblIds which I am trying to convert to gene symbols using BioMart. Some of the ensemblIds will not be mapped, so I need to subset the rows of the matrix. It seems like BioMart adds one ensemblId or at least I see one more result ultimately than should be. Here is what I am trying to do:
library(dplyr)
library(tibble)
ids <- rownames(mtx_combined)
mart <- useDataset("mmusculus_gene_ensembl", useMart("ensembl"))
list <- getBM(filters= "ensembl_gene_id", attributes= c("ensembl_gene_id",
                           "mgi_symbol","description"),values=ids,mart= mart)
mtx_mapped <- as.matrix(as.data.frame(mtx_combined) %>%
rownames_to_column("gene") %>%
filter(rownames(mtx_combined) %in% list$ensembl_gene_id) %>%
column_to_rownames("gene"))

rownames(mtx_mapped) <- list$mgi_symbol   //Fails

Length of the returned list is 18337, number of rows in the initial matrix is 18376. After subsetting I expect to get 18337 rows, but getting one less - 18336. What am I doing wrong here?

What is the extra line? Just print them out and see. Could it be adding a header, for example?

There’s no 1:1 mapping between the IDs of different database schemas. As a consequence, some Ensembl gene IDs map to multiple MGI symbols, or to none (and vice versa).
Therefore, you can’t assume that the number of result rows will be equal to the number of query IDs.


text:
If I modify a PDB file with a specific mutation, how to minimize energy?
I have a series of proteins that are all phylogenetically related. Only one of this proteins is currently X-ray crystallized.
Is it valid to modify the reference PDB file (I mean by hand) to match the mutated sequence?
And if so, is there any way to minimize the protein energy to find the most likely structure? (Trying to match as closely as possible the reference structure)

I would recommend you try the FoldX mutation engine.  It is very easy to use and state of the art (it is the main competitor to expensive commercial software).  It can do both of what you want - the mutation and energy optimisation.  It may not work as well as molecular dynamics software, but for your simple task, it is a reasonable solution.
A simple command line can do it:
First, repair your structure (since, usually PDBs are not energetically optimised, so you could get better energy anyway without doing this):
FoldX --command=repairpdb --pdb=1AJY.pdb

Then, either use AlaScan or PositionScan to mutate.  AlaScan only mutates to Alanine.  PositionScan can mutate to any amino acid. More advanced mutations can be done with the BuildModel command.
FoldX --command=PositionScan --pdb=1AJY_Repair.pdb --output-file=1AJY_Repair_MA30d --out-pdb=false --positions=MA30d

Note: 'd' means try mutate the residue to 24 different amino acids; or, you can specify any of those 24, including the standard 20.
Which produces two files:

PS_1AJY_Repair_MA30d_scanning_output.txt
energies_30_1AJY_Repair

(It can also output the mutated pdb file by setting --out-pdb=true.)
Either of those files can be used to determine change in free energy upon mutation (ddg ~ delta delta g).  The second file has additional information (the individual values used to calculate total change in free energy).
Note: The FoldX website has very good documentation for each FoldX command.  You can find any additional information you require there.

Yes you can modify the reference PDB file and look for the changes and for this purpose you need visualizers. One among these is Chimera. You can easily carry out the energy minimization steps using Chimera, first mutate your PDB file manually and then provide it as an input for Chimera and the steps for doing energy minimization are fairly simple as given in Chimera docs here.

The previous answers suggest Chimera and FoldX, which require a license for commercial use.
Commercial users can use pdbfixer or a similar tool to perform mutations (using e.g. applyMutation in the API).
You can then use OpenMM to refine the mutated structure.
Granted, this requires some comfort writing Python scripts to use the APIs.


text:
differentially expressed intron analysis?
I have used Deseq for estimating the differentially expressed introns for the RNA seq data using the UCSC intron bed file. can Deseq be used for differentially expressed intron analysis?

Welcome to the site. You write Deseq in the question but you tagged with deseq2 (of the DEseq2 package). Which one do you really use? Also if you did it can be used, but it seems like you have problems with the results or something, if this is the case, you will need to provide more of information like the relevant part of the code used

Typical tools to use for this would be the DEXseq packages (as done here) or rMATS, which will calculate intron retention rates.
I would generally suggest not summing intronic counts by gene and using that in DESeq2 (I assume that's what you meant by DESeq, since there's no reason to use the original DESeq for anything). What you'd be getting if you did this is less expressed introns and more sample-specific bias in enriching for mature versus pre-mRNA.

Thank you Ryan. Yeah, I will use either DEXseq or rMATS in that case.


text:
tximport txidcol and geneidcol selecting wrong column
I'm trying to run a differential expression analysis of 2 samples of the same species, but with different treatments using DESeq2. I'm using Ensembldb to create a tx2gene dataframe for tximport, however tximport does not select the columns I specify, and instead always selects the tx_biotype leading to the following output:
.                    25_quant    28_quant
nontranslating_CDS    10422.41    10022.74
protein_coding     22196722.59 24002269.26

I'm running tximport using the following params:
txi_tx <- tximport(files, type = "salmon", tx2gene = tx2gene,
                   txIdCol = "tx_id", geneIdCol = "gene_id")

I've gotten around this by deleting all columns from my tx2gene DF apart from tx_id and gene_id. However, I think this might be causing problems in the DESeq2 analysis somehow, and I would like to rule out that possibility.
This is what my names(tx2gene) looks like:
"tx_id"            "tx_biotype"       "tx_seq_start"     "tx_seq_end"
"tx_cds_seq_start" "tx_cds_seq_end"
"gene_id"          "tx_name"

What am I doing wrong?

By setting type="salmon" you're overriding your settings for things like txIdCol. The various type options (except none) are presets to set things correctly for the output from various programs like salmon or kallisto. If you need to specify custom column names, set type="none".
As an aside, if I'm interpreting your question correctly you have only two samples in total. You can not perform a meaningful analysis with an unreplicated experiment, so if that is indeed your case then don't waste your time further with this dataset.

Thank you, that explains it. As for the lack of replicates this is for an undergraduate research project, and we're just trying to get our R code to work until we get actual data.

Ah download a public dataset and try to reproduce a published analysis with that. If you can set similarish (or better!) results as those published then you're good to go :)


text:
Behavior of `--reference` flag with samtools sort
I'm about to run samtools sort (version 1.6) and I'm reviewing the available configuration options.
$ samtools sort
Usage: samtools sort [options...] [in.bam]
Options:
  -l INT     Set compression level, from 0 (uncompressed) to 9 (best)
  -m INT     Set maximum memory per thread; suffix K/M/G recognized [768M]
  -n         Sort by read name
  -t TAG     Sort by value of TAG. Uses position as secondary index (or read name if -n is set)
  -o FILE    Write final output to FILE rather than standard output
  -T PREFIX  Write temporary files to PREFIX.nnnn.bam
      --input-fmt-option OPT[=VAL]
               Specify a single input file format option in the form
               of OPTION or OPTION=VALUE
  -O, --output-fmt FORMAT[,OPT[=VAL]]...
               Specify output format (SAM, BAM, CRAM)
      --output-fmt-option OPT[=VAL]
               Specify a single output file format option in the form
               of OPTION or OPTION=VALUE
      --reference FILE
               Reference sequence FASTA FILE [null]
  -@, --threads INT
               Number of additional threads to use [0]

I'm curious about the --reference flag. I've never used it before with samtools sort and I'm wondering what behavior it invokes. The online documentation doesn't even mention the flag. Under what circumstances would I want/need to specify the reference genome?

You would only need to supply a reference when using CRAM output encoding. The CRAM format allows storing entries without explicitly encoding the base calls, since one can instead encode the differences to the reference genome.
With BAM output the --reference option has no effect.


text:
GATK CombineVariants complains the contig order in the VCF files
I have called variants on two strains of C. elegans separately. I now want to merge the VCF files into one using the following code:

Create a sequence dictionary of the reference sequence
Sort the VCF files with Picard
Merge the sorted VCF files using GATK

picard CreateSequenceDictionary \
    REFERENCE=c_elegans.PRJNA13758.WS263.genomic.fa \
    OUTPUT=c_elegans.PRJNA13758.WS263.genomic.dict

picard SortVcf INPUT=strain1.vcf \
    OUTPUT=strain1sorted.vcf \
    SEQUENCE_DICTIONARY=c_elegans.PRJNA13758.WS263.genomic.dict

picard SortVcf INPUT=strain2.vcf \
    OUTPUT=strain2sorted.vcf  \
    SEQUENCE_DICTIONARY=c_elegans.PRJNA13758.WS263.genomic.dict

GATK --analysis_type CombineVariants \
    -R c_elegans.PRJNA13758.WS263.genomic.fa \
    --variant strain1sorted.vcf \
    --variant strain2sorted.vcf \
    -o all.vcf \
    -genotypeMergeOptions UNIQUIFY

The last command gives me the following error message:
ERROR MESSAGE: Input files variant and reference have incompatible contigs. Please see https://www.broadinstitute.org/gatk/guide/article?id=63 for more information. Error details: The contig order in variant and reference is not the same; to fix this please see: (https://www.broadinstitute.org/gatk/guide/article?id=1328),  which describes reordering contigs in BAM and VCF files..
##### ERROR   variant contigs = [I, II, III, IV, MtDNA, V, X]
##### ERROR   reference contigs = [I, II, III, IV, V, X, MtDNA]

But I have sorted the VCF files using Picard, so I don't know what else to do.
Your help is appreciated.
Update: I asked this question on the GATK forum on 4th May 2018 (today).

GATK cannot combine VCF files generated by Sniffles, the variant caller that I used to call structural variants (answer from the GATK team).

Answered in the GATK forum

Try regenerating the index files for your VCFs. Picard SortVcf doesn't do it for you iirc, so when GATK looks at the index files (before opening the VCFs themselves) it sees the old order and complains.

Hi, thanks for the answer, it's already helpful. However, It would be even nicer for community to write the actual answer and to provide the link just as the source (you can [edit] your answer).


text:
Create gene network from list of genes
This has also been posted on biostars
I had an initial bipartite transcriptional regulatory network for bees. My idea is to create a similar network for termites so as to check if similar genes are hubs in both the networks (being that both are eusocial). For this purpose, I took the sequences of the genes in the honey bee GRN and did a BLAST to identify similar genes in the species of eusocial termite that I am interested in. The problem is that I am out of ideas on how to infer a network from the set of genes that I have.
I cannot use Genemania which can only be used for genes from specific species like humans, flies etc. I cannot use STRING because that uses protein sequences and the species I work with is also not on it. Most tools I found can be used for network creation only if the genes I have is from well studies species like mus musculus, drosophila etc.
Is there anyway in which I can create a network from the set of genes that I have?
I do not have interaction information. I have a network of genes from honey bees. I did a BLAST of the genes in that network and retrieved the gene names for the related genes that I got from a species of termites. I would like to create a gene network of that to check if the same genes form the hubs in both the networks.
The gene network in bees is a bipartite transcription regulatory network. I haven't checked the gene ontology of the termite genes but I would think that the genes should be having similar functions because they are BLAST results. There are around 2200 genes in the honey bee GRN and I could find approx 1000 genes from the BLAST and UNIPROT for termites. The honey bee GRN was generated experimentally by a different lab so I could not replicate what they did with ARACNE.
What do you suggest that I do to find similar genes then if a BLASTN is not sufficient information?


text:
BLAST(n): No hits found
I am currently exploring the BLAST program, just for testing purpose i generated two FASTA files, containing two genes A and B, such that B is just a motif of repeated 'G's that occurs in A.
file A.fna:
>A
ACCAACAGACGGGGGGGGGGGGGAGAAAACCGGCCCCCGCAGGGACGCCACCCAGAGGGAAGGG

file B.fna:
>B
GGGGGGGGGGGG

So by my understanding BLAST should find the local alignment of B in A easily? But the program returns 'No hits found'? Maybe I am doing something wrong here?
I executed the following:
makeblastdb -in A.fna -dbtype 'nucl'
makeblastdb -in B.fna -dbtype 'nucl'
blastn -db A.fna -query B.fna

output of blastn:
BLASTN 2.7.1+

Reference: Zheng Zhang, Scott Schwartz, Lukas Wagner, and Webb Miller (2000), "A greedy algorithm for aligning DNA sequences", J Comput Biol 2000; 7(1-2):203-14.

Database: A.fna
           1 sequences; 64 total letters

Query= B

Length=12

***** No hits found *****

Lambda      K        H
    1.33    0.621     1.12

Gapped Lambda      K        H
    1.28    0.460    0.850

Effective search space used: 480

  Database: A.fna
    Posted date:  May 3, 2018  5:17 PM   Number of letters in database: 64   Number of sequences in database:  1

Matrix: blastn matrix 1 -2 Gap Penalties: Existence: 0, Extension: 2.5

Thanks for any suggestions (:

I'm not sure but could it be that there are too many gaps and then it doesn't match it. BLAST doesn't just search for a sequence inside another one, but looks for similarities between sequences, and these two sequences are not similar (even if one contains the other)

But is that not exactly what BLAST is supposed to do (local alignment)? There is a perfect local alignment (without any gaps). In terms of global alignment the two sequences are very different (and gap penalties are huge), but here we are talking about local alignments !?

Yes, it should, glad we have terdon to answer the question

There are three possible problems that come to mind.

Masking
Blast will mask low complexity regions by default. Since your sequence is nothing but Gs, it is a safe bet that it is being masked, so no hits will be found for it.

Score/e-value thresholds
Another source of complication is that even if a match is found, that match will have very bad scores. Both the actual score of the alignment and the e-value will be very bad. Since this is such a simple sequence, it will always score badly.

Word-size
The way blast works is (simplifying a bit) by finding a match for N residues (the word size) and then extending that match if extending increases the score. If your query sequence is shorter than the word size, no match will ever be found.

With all that said, the specific problem in your case seems to be a combination of the word-size and masking. If you run your query setting the word-size to 11 and disabling masking, as shown below, you should get results.
blastn -db A.fna -query B.fna -word_size 11 -dust no

That should return two hits (two because the word size is one less than the query size, so it can align perfectly at two locations):
> A
Length=64

 Score = 23.3 bits (12),  Expect = 5e-05
 Identities = 12/12 (100%), Gaps = 0/12 (0%)
 Strand=Plus/Plus

Query  1   GGGGGGGGGGGG  12
           ||||||||||||
Sbjct  11  GGGGGGGGGGGG  22

 Score = 23.3 bits (12),  Expect = 5e-05
 Identities = 12/12 (100%), Gaps = 0/12 (0%)
 Strand=Plus/Plus

Query  1   GGGGGGGGGGGG  12
           ||||||||||||
Sbjct  12  GGGGGGGGGGGG  23

Note that the default word size for BLASTn searches is usually 11. I don't know why the blast+ version (2.2.31+) I have on my Ubuntu machine didn't use that as a default. The manual says it should. But I got no hits when I didn't explicitly set the word size to 11.

that makes a lot of sense (i was not familiar with the concept/term of masking) thanks (:

I came across this thread while searching for a solution to a similar problem I was experiencing. I am inputting markers that are 20 basepairs in length as the query sequence and haven't been getting hits. As terdon mentioned in their excellent reply, setting word size to something smaller like 11 helps for these shorter query sequences.
I didn't get any results on my 20 bp query sequences when using blastn against the sunflower genome, but I got plenty when setting word size to 11. So, why does setting word size for blastn to 11 improve the results despite documentation stating that default word size is 11? Well it turns out that this is very misleading. This webpage explains that the default algorithm used in blastn is actually megablast, not blastn. The default word length for this default megablast algorithm is actually 28!
The documentation is quite unclear about this and many websites providing tutorials for blastn also seem to miss this key detail. Hope this helps anyone else who struggles with this!

Your hyperlink doesn't seem to be working any longer, but megablast vs blastn is a huge difference. Command line -task blastn or python task="blastn" will switch the algorithm.


text:
Subset on multiple genes in Seurat
I know that I can do subsetting on just one gene in Seurat:
seurat_subset <- SubsetData(seurat_object, subset.name = neuron_ids[1], accept.low = 0.1)

However, I want to subset on multiple genes. Is there a way to do that? I just do not want to do manual subsetting on 10 genes, then manually getting @data matrix from each subset, and recreating seurat object afterwards.

I was able to achieve this in the following way:
require(data.table)
cell_names <- vector(mode="character")
for (i in 1:length(neuron_ids)) {
  seurat_subset <- SubsetData(seurat_object, subset.name = neuron_ids[i], accept.low = 0.1)
  genes <- colnames(seurat_subset@data)
  cell_names <- c(cell_names, genes[! genes %chin% cell_names])
}
seurat_subset <- SubsetData(seurat_object, cells.use = cell_names)

Would be interesting to know if Seurat provides such functionality out of the box.

In Seurat V3 you should use seurat_subset <- SubsetData(seurat_object, cells = cell_names) as cells.use is no longer supported.


text:
Error filtering SNPs with Plink
Im trying to filter SNPs using PLINK with this command. Also, the .map and .ped files were converted to a plink format from a vcf file (python executable.py --input allindividuals.vcf)
plink --ped allindividuals.ped --map allindividuals.map --allow-extra-chr --indep-pairwise 50 10 0.1

but at the time I count the SNPs in the file plink.prune.in it has only 15 SNPs and in the file plink.prune.out it remains millions of them. In addition, when I change the parameters, the result is the same, no matter how many times I change it.

Answered this on Biostars: the issue was that there were only 14 samples, and almost no genetic variation between them.

@ChristopherChang if you have the time, it would be great if you could post your answer  here too!

@terdon here is the link: https://www.biostars.org/p/281334/#313252


text:
Does the galGal5 chicken assembly have a chromosome 29?
The chromosome sizes at UCSC don't seem to contain chr29:
ftp://hgdownload.soe.ucsc.edu/goldenPath/galGal5/bigZips/galGal5.chrom.sizes
It has a chr28 and a chr30. Am I missing something or is there some piece of history that led to the omission of this chromosome?

The current ensembl entry doesn't have a 29 either.  The archived ensembl assembly  lacks 29 30, and 31 and 33 and LGE64.
The chromosomes after 30 are tiny, so they might not be visible in a karyotype.  They probably realized that "chr 29" was really attached to some other chromosome.


text:
BioMart Ensembl id imperfect retrieval
I have a list of gene symbols:
c("cd45", "Tmem119", "CD11b", "P2Yr12", "Csf1r", "Bst2", "Cd74",
"Cx3cr1", "Trem2", "Lyz2", "GLAST", "GFAP", "ALDH1L1", "Aldoc", "Aqp4",
"Glul", "S100a", "Olig1", "Olig2", "Olig3", "Mbp", "Pdgfra", "Pecam",
"Cldn5", "Cldn10", "Epas1", "Crip1")

If I feed it to BioMart to get ensembl ids:
mart <- useDataset("mmusculus_gene_ensembl", useMart("ensembl"))
list <- getBM(filters= "mgi_symbol", attributes= c("ensembl_gene_id",
                           "mgi_symbol","description"),values=symbols,mart= mart)

the following genes are missed:
"cd45"   "cd11b"  "p2yr12" "glast"  "s100a"  "pecam"

All of them are pretty well-known genes, and I can manually find their ensembl ids by googling it, for instance:
http://www.informatics.jax.org/marker/MGI:97810
I tried supplying aliases, but it does not change the output. So, from my understanding BioMart is not working properly either because I am doing something wrong here or because BioMart itself is not a good tool to use. Is there a better way of getting the mapping that would map all of the gene symbols?

The genes that are missed are probably not official mgi symbols.
You might wanna look them up at mgi:
cd45 -> Ptprc

cd11b -> Itgam

p2ry12 -> P2ry12 (?) (CAPS sensitive?)

glast -> Slc1a3

s100a -> S100a1

pecam -> Pecam1 (?)

I suspect biomaRt is not the problem here, the only better way I can come up with is to download the ENSEMBL annotation (GTF) file, and get the symbol out of this annotation with AWK or in R.
Here an example in R, since you don't want to use AWK:
file <- read.table("Mus_musculus.GRCm38.78.gtf", sep = "\t", comment.char = "#")

exons <- subset(file, V3 == 'exon')

V9 <- file[,"V9"]

ensembl_names <- sub("gene_id ","",V9)

split <- strsplit(ensembl_names, ";")

splitting <- sapply(split, function(x) x[1])

split_symbol <- strsplit(ensembl_names, "gene_name ")

splitting_symbol <- sapply(split_symbol, function(x) x[2])

splitting_symbol2 <- strsplit(splitting_symbol, "; ")

symbol_names <- sapply(splitting_symbol2, function(x) x[1])

ensemble_symbol_big <- cbind(splitting, symbol_names)

ensembl_symbol <- unique(ensemble_symbol_big)

colnames(ensembl_symbol) <- c("ENSEMBL", "Symbol")

write.table(ensembl_symbol, "ensembl_symbol.txt", sep = "\t", row.names = F)

I am not sure what you mean by look up at mgi. I tried supplying mgi instead of mgi_symbol, and it did not work. Not sure whether it is an issue with biomart or not. How to do that with AWK? Can you give any links for doing that, tutorials? To learn whole new language just to do translation of ensembl ids to symbols seems like an overkill a bit, if it is not very easy.

What I mean with "look up at mgi", is look up the genes that are missed: e.g., cd45, and cd11b. You'll see that your names do not match the name given after Symbol, so how do you expect biomaRt to recognize them? I think AWK is a great tool for bioinformaticians, but if you don't want to make use of this tool, don't. I do not want to force you into something.


text:
Average of logCPM
I am a bit in doubt how to calculate correctly the mean (average) in the logarithmic scale. I know e.g., that a division in the normal scale becomes a subtraction in the logarithmic scale, and so forth, that's why I was wondering what the correct way is to calculate the mean (since the mean is calculated by additions and divisions in the normal scale).
So for example, if I have in RNAseq data a group with 3 samples, and have log2 CPM values for each sample. How can I calculate the mean of the group correctly?
Do I first need to get CPM (without log2), calculate the mean and then apply log2? Or can I just ignore the fact that the log2 CPM values are in the log scale, and just calculate the mean as if it was the normal scale (add all up and divide by n)?

You can calculate the mean of a value that is a log2 of another value as you would normally. However mean(log2(CPM), log2(CPM), ...) != log2(mean(CPM, )), choosing one or another depends on what do you want to do.

And which one to choose if you want to represent the average of the group?

I usually use mean(log2(, as the log2CPM are the normalized values, so the "important" ones.

The mean of log-transformed values is the log-transformed geometric mean of the untransformed values (i.e., $log_{2}(geometric~mean)$). This can be useful when the values you want to summarize have vastly different ranges (so one doesn't dominate the average). I would expect that to generally NOT be the case for CPMs of the same gene, so it's unclear to me what the point of that would be. It can also be difficult to interpret these results due to the geometric mean being dependent upon the variance:
> originalValues = rnorm(1000, mean=10, sd=c(rep(1, 500), rep(2, 500)))
> mean(originalValues[1:500])
[1] 9.991464
> mean(originalValues[501:1000])
[1] 10.04407
> log2(mean(originalValues[1:500]))
[1] 3.320696
> log2(mean(originalValues[501:1000]))
[1] 3.328272
> mean(log2(originalValues[1:500]))
[1] 3.313263
> mean(log2(originalValues[501:1000]))
[1] 3.300103

Above I created two groups of 500 values each. Each group has approximately the same mean (10) but with different variance (1 vs 2). As you can see, the second group of 500 values randomly has a slightly higher average. The log2(mean) preserves this difference, but the log2(geometric mean) inverts the relationship such that group 2 now appears to have a lower mean value. This isn't really wrong, it's the difference between mean and geometric mean, but it's something that needs to be kept in mind.
This has implications if your want to compare values between genes, since your comparison metric is then not just dependent upon the mean, but also the variability of the data (such that more variable genes will have systematically lower geometric means).
In short, the geometric mean is useful when values are on different scales (such as CPM values of all genes in a sample, but probably not CPM values of the same gene between samples). If that's not the case, then you're mostly removing interpretability from the results.


text:
How to enlarge a section of y-axis?
I am using ggplot2 to plot a scatter plot.
library(ggforce) # required by facet_zoom

g <- ggplot(data.plot, aes(x = Methyl_Average, y = Average_Gene_Expression))
g1 <-  g + geom_point(color = "steelblue") + theme_bw(base_family = "Times") +
       ggtitle( "Scatter Plot") +
       theme_bw(base_family = "Times") +
       theme(plot.title = element_text(hjust = 0.5))

As most of my data points are in the range of 0 - 200 on y-axis, I want to zoom this section. I tried doing it using facet_zoom from ggforce package :
g1 + facet_zoom(y = Average_Gene_Expression > 0 & Average_Gene_Expression < 200)

Instead of two plots in one figure, I want one plot (like first figure), where 80% of y-axis displays the data points in the range of 0 - 200 and 20% of the y-axis displays the remaining and x-axis unaltered. Which function/package should I use do it in R?

I'm voting to close this question as off-topic because it is not about bioinformatics, it is about representing results. Which might be better suited in academic.stackexchange.com

I voted to leave it open because I think visualization is an integral part of bioinformatics and the presented data to be visualized are obviously very domain-bioinformatic specific.

As per this link this is not possible with ggplot2. Also check this answer.

I would strongly discourage you from making discontinuous axis, it's going to be very confusing for a reader.
The facet plot you proposed seems like a good solution to me. Alternatively you can use log transformation. To demonstrate I made it on simulated data that look appox like yours :
set.seed(940401)

data.plot <- data.frame(Methyl_Average = c(0.05, rbeta(699, 2, 5) + 0.15, runif(300, 0.18, 0.85)),
                        Average_Gene_Expression = c(165, exp(rnorm(999, 3))))

library(ggplot2)
require(gridExtra)

g <- ggplot(data.plot, aes(x = Methyl_Average, y = Average_Gene_Expression))
g1 <-  g + geom_point(color = "steelblue") + theme_bw(base_family = "Times") +
       ggtitle( "Scatter Plot") +
       theme_bw(base_family = "Times") +
       theme(plot.title = element_text(hjust = 0.5))

sc_plot <- g1 +
           labs(x = "Differential Methylation", y = "Differential Expression")
logsc_plot <- g1 +
              labs(x = "Differential Methylation", y = "Differential Expression (log10 scale)") +
              scale_y_continuous(trans='log10')

grid.arrange(sc_plot, logsc_plot, ncol = 2)

Nice answer. One minor point: the numbers on the y-axis of the right plot is still DE, not log10(DE). The scale is log transformed but the numbers aren't. I'd call it Differential Expression (log10 scale) or something.

fair point, but I am too lazy to redo the image, I changed it in the code...


text:
How to get rid of the temp files while using "featureCounts" for extracting readcounts from bam files?
I used this command for extracting read counts for almost 60 samples using featureCounts.
featureCounts -a /gencode.v27.primary_assembly.annotation_nochr.gtf -F GTF -p -s 2 -T 8 -o counts.txt sample.sorted.bam

I see many "temp-core-100373-A4BF0116E148.sam-TH00-BK004414.tmp" files like this dumped in my directory. I'm able to delete them only manually. Is there a way to get rid of these tmp files from my directory while using the command?

You shouldn't delete them while featureCounts is running, it's creating read-name-sorted temporary files so it can more easily look at both mates together. Just wait for it to finish running and it will clean up the files itself. Doing so beforehand is liable to cause an error.

Thanks. But will there be any memory problem for those temp files? Asking this cz right now I have only 60 samples, for those samples I removed temp files manually. It didn't clean up even after finishing the job. I have to use featureCounts on at least 900 samples in future.

There aren't memory issues involved, though I suppose you could run into temporary space constraints. If featureCounts didn't clean up after itself then I encourage you to check and ensure it didn't encounter an error. I have never had it leave temp files around otherwise.


text:
Reads mapped to exonic, intronic and intergenic regions
After the alignment step I checked the rnaseq metrics of all the samples. Among 40 samples three samples show high percentage of reads mapped to intronic regions. What could be the reason?
Samples   Exonic            Intronic    Intergenic
Sample1 545479 (12.8%)  2512309 (58.8%) 1217201 (28.5%)
Sample2 372836 (8.3%)   2556934 (56.8%) 1573032 (34%)
Sample3 529618 (7.3%)   3934615 (54.5%) 2750720 (38.1%)

I also see that 35-40k genes among 58k genes having "zero" read counts. Is this due to contamination? What could be the reason? What does reads mapping to exonic, intronic, intergenic tell?
Update: I used hisat2 for alignment. I used human genome "grch38_snp_tran" from hisat2 website. Libraries were generated using ribosome depletion kit.

How did the overall alignment rate of those samples compare to the others that had more reasonable intronic alignment rates?

I see that all the 37 samples have overall alignment rate between 90-96% and for those 3 samples it is 69%, 79%, 70% overall alignment rate.

Sounds like you have degraded samples

U mean those 3 samples are degraded? What I should do now? Should I remove them for further analysis?

It’s a likely cause, but regardless you should exclude them.

Sure. the high percentage of mapping to intronic region is due to degradation?

Could you add the quality metrics of these files and other of the same run to the question?

There is the potential that these samples had low quality RNA or genomic DNA contamination (depending on the sample type and RNA isolation method)

What software did you get this output from?

As @DevonRyan mentioned, it's very likely that those samples were degraded, which is good justification for excluding them from subsequent analysis.

I would expect at least 30% of reads from a total-cell, ribo-depleted RNA-seq to be exonic. Less suggests something when wrong.
As well as degradation, another explanation would be contamination with genomic DNA.
RNA obtained from the nuclear or cytoplasmic fraction of the cell might have a exonic content lower than 30%, while poly-A selected RNA would be expected to have higher, more like 66%.


text:
How to decide number of neighbors and resolution for Louvain clustering
I am using Louvain clustering (1,2) to cluster cells in scRNAseq data, as implemented by scanpy.
One of the parameter required for this kind of clustering is the number of neighbors used to construct the neighborhood graph of cells (docs).
Larger values result in a more global view of the manifold, leading to lower number of clusters, while reducing the number of neighbors goes in the opposite direction. However, it is unclear how to choose this parameter.
The resolution parameter seems to work in the opposite way.
Do you know of any methodology and/or rule-of-thumb to define these parameters? E.g. depending on the size of the dataset?

Levine, Jacob H., et al. "Data-driven phenotypic dissection of AML reveals progenitor-like cells that correlate with prognosis." Cell 162.1 (2015): 184-197.
Blondel, Vincent D., et al. "Fast unfolding of communities in large networks." Journal of statistical mechanics: theory and experiment 2008.10 (2008): P10008.

The Louvian algorithm has an issue with disconnected communities. The Leiden algorithm is recommended to resolve this. See Traag et al. https://arxiv.org/abs/1810.08473

If I remember correctly the phenograph paper shows that their algorithm is fairly robust for different (reasonable) values for k.

@gc5 did you get the answer? I also want to know what is the best resolution ? As we can get more communities when we increase resolution and vice-versa. Thanks

@KhalidUsman no satisfying answer so far.

A general rule of thumb is that in order to improve the variance $n$ times you need $n^2$ neighbours. This is only applicable if you consider the $n^2$ nearest neighbours of a cell to be biologically identical (i.e. "similar enough"); if your data includes 10 types of cells with 10 cells each, then using the 20 nearest neighbours for smoothing will obscure the data.
As far as I know, there is no single best answer to this question. I would suggest trying different numbers and sticking to what agrees more with the biology of the dataset.

Can you put some references on the rule of thumb you wrote? However, I ended up selecting progressive resolutions. Of note, different parts of a dataset may need different resolutions. Thanks.

I read a paper where the authors used different combinations of parameters to iteratively perform LSI until some criteria. Maybe a similar strategy can be applied to Louvain clustering as well?
The paper can be found here: https://www.nature.com/articles/s41587-019-0332-7.
You can also check their code to see how they implemented it.


text:
Manually define clusters in Seurat and determine marker genes
I want to define two clusters of cells in my dataset and find marker genes that are specific to one and the other. Is there a way to do this in Seurat? Say, if I produce two subsets by the SubsetData function, is there a way to feed them into some other function that would calculate marker genes? If not, what other packages would you recommend for doing that?
If you look here:
https://satijalab.org/seurat/seurat_clustering_tutorial_part2.html
I just need a way to define ident myself, the number of levels (2) and assign numbers to each cell (0, 1), and then run DE between 0 and 1 clusters which is obvious how to do afterwards.

A related issue was answered on github. Just as an alternative reference to the excellent answers below)

I think you are looking to FindAllMarkers function from Seurat. As you said, you just have to define your ident, that have to have the structure of a table (cell names as names and cluster as value):
pident=as.factor(clusters)
names(pident)=cellNames
object1@ident=pident

And then run the FindAllMarkers function:
FindAllMarkers(object1, min.pct = 0.25, min.diff.pct = 0.25)

You can specify several parameters in this function (type of DE to perform, thresholds of expression, etc).

Seurat has functions for adding metadata and setting identities. Get unique cell names:
cell.labels <- seuratobject@ident

Replace column and its name with your cluster labels (e.g.), then:
seuratobject <- AddMetaData(seuratobject, metadata=cell.labels)
seuratobject <- SetAllIdent(seuratobject, id='yourclusterlabels')

Because you want to contrast two clusters against each other, I suggest using FindMarkers() as opposed to FindAllMarkers():
FindMarkers(object, ident.1, ident.2)

It can also compare combinations of clusters.


text:
Transform and feed data into DESeq2 with DESeqDataSetFromMatrix
There is a normalized expression matrix. I split it into two and want to do DE on the two cells' subsets. I am having trouble transforming it into the format that DESeq2 would accept. The function that I would think I need to use is the following:
dds <- DESeqDataSetFromMatrix(countData = cts,
                          colData = coldata,
                          design= ~ batch + condition)

It would be perfect if I could somehow feed two matrixes separately, but I can not find it. Any suggestions would be greatly appreciated.

Update

I tried the following:
cells <- c(cell_names, cell_names_cyt)
groups <- c(rep("NE", times=length(cell_names)), rep("EX",
times=length(cell_names_cyt)))

colData <- data.frame(
  Cells <- cells,
  Groups <- groups
)

dds <- DESeqDataSetFromMatrix(countData = df,
                          colData = colData,
                          design= ~ Groups)

But getting the error:

Error in DESeqDataSet(se, design = design, ignoreRank) :
    some values in assay are not integers

Of course, they are floats since I already normalized the matrices.

Update

I did the following since I am using Seurat:
df <- as.data.frame(as.matrix(seurat_object@raw.data))
df <- df[,colnames(df) %in% cells]

But then the following:
dds <- DESeqDataSetFromMatrix(countData = df,
                          colData = colData,
                          design= ~ Groups)

Gave me an error:

Error in DESeqDataSet(se, design = design, ignoreRank) : all variables in design formula must be columns in colData

Don't input normalized counts, use raw counts.

Yeah, I used raw counts retrieved from Seurat, but then got another error.

Error in DESeqDataSet(se, design = design, ignoreRank) : all variables
  in design formula must be columns in colData

I think that means just what it says...you don't have a "Groups" column in coldata.  Did you spotcheck your ColData data frame to make sure it looks the way it should?
> cells <- c("1","2","3")
> groups <- c("a", "b", "c")
> test <- data.frame(Cells <- cells, Groups <- groups)
> test
  Cells....cells Groups....groups
1              1                a
2              2                b
3              3                c
> colnames(test)
[1] "Cells....cells"   "Groups....groups"


text:
What kind of analysis can be done with differential expression of transcription factors?
I have two different stem cell types and their respective gene expressions from RNAseq.
I noticed that there is differential expression in some of the genes coding for histones as well as transcription factors?
What kind of interesting downstream analysis might be available?
Is it feasible to survey ENCODE for downstream TF binding sites?
Can anything interesting be investigated with regard to the histone gene data?
Thanks in advance.

TF usually bind in the upper region of the genes it regulates, so you could see which genes are in the region below the binding site. For histone related data, I'm not aware of how it work but perhaps you can find some leads looking at methylation studies.

Expanding my comment.
TF usually bind in the upper region of the genes it regulates, so you could see which genes are in the region below the binding site. Perhaps there are some databases for TF.
I found a website base on the ENCODE project with the genes regulated by each TF.You can use them for pathway or gene set enrichment analysis, but I am not aware of any specific analysis for genes regulated by a single TF.
The histone data is harder if it was methylation I would say to look if different states of the histones are in different spatial state, but with gene expression is harder to conclude anything. Nevertheless you can check methylation studies to see if gene expression is associated with methylation changes.

It sounds like you’ve already done differential expression analysis (with genes or transcripts) and a pathway enrichment analysis to find these. Differential expression of these functional groups is an interesting finding in its own right but there are several further investigations that could follow, depending on your biological research interests.
Clustering the expression patterns of the differentially expressed genes would be able to show if there are relationships between them (such as inter-gene correlations) or relationships to other phenotype data. This is commonly performed with unsupervised hierarchical clustering and displayed by a heatmap. You could also look into whether these differentially expressed genes are expressed in a tissue-specific manner. There is plenty of reference material available from ENCODE, FANTOM4/5, or GTEx that could be used for comparison (as well as cancer datasets) and the Human Cell Atlas Project (currently underway).
Another approach is pathway structure. There are many hierarchical annotation of gene functions such as gene ontology (GO terms) that overlap and are ranked: identifying smaller more specific pathways for future investigations. No pathway annotation is perfect so it’s best to replicate important findings across pathway databases (ideally different expression datasets as well). These functional groups of genes aren’t unordered: the structural relationships between pathways are well characterised in some cases. This includes many gene regulatory networks involving transcription factors. These graph structures are supported on databases such as WikiPathways and Reactome. Network analysis is relatively unexplored compared to differential expression.
While this is a bioinformatics question, there’s always a need to consider experimental validation as a future direction. If you already have a hypothesis based on your analysis, it may be more appropriate to design an experiment or contact experimental collaborators. If you have access to further genomics data, you could also explore how these differences are regulated: CAGE-Seq for enhancers and lncRNA, ATAC-Seq, and 3C/4C for chromatin structure, or Bisulfite-Seq and ChIP-Seq for epigenetic marks. An integrative analysis of expression and epigenomics data would be especially compelling and would likely account for different activity of transcription factor binding sites. Unfortunately, this can’t be done without more ‘omics data available but it would greatly benefit investigations of histones if it is feasible.


text:
get gene lines from gtf file
I would like to retrieve gene lines from a GTF file for which I only have exons & transcripts lines (output from Cufflinks) and alternative splicing possible. I need gene lines for compatibility with a pipeline dealing with GTF in Ensembl format.
Ideally, I would have gene lines representing the longest transcript - but I am open to discussion about that. The aim would be to separate the genome in genic & intergenic portions - boundaries being defined by start-end coordinates of the genes.
Sample input (e.g. gene "CUFF.105" with 5 transcripts):
cat subset_105.test | awk '$3=="transcript"{print $0}'

tig00000046 Cufflinks   transcript  26170   42766   202 +   .   gene_id "CUFF.105"; transcript_id "CUFF.105.3"; FPKM "0.2320304094"; frac "0.197788"; conf_lo "0.164885"; conf_hi "0.299176"; cov "5.188855";
tig00000046 Cufflinks   transcript  26170   42766   266 +   .   gene_id "CUFF.105"; transcript_id "CUFF.105.1"; FPKM "0.3061260755"; frac "0.249779"; conf_lo "0.228860"; conf_hi "0.383392"; cov "6.845843";
tig00000046 Cufflinks   transcript  26170   39469   470 +   .   gene_id "CUFF.105"; transcript_id "CUFF.105.2"; FPKM "0.5403628578"; frac "0.161731"; conf_lo "0.372512"; conf_hi "0.708214"; cov "12.084038";
tig00000046 Cufflinks   transcript  28928   39469   1000    +   .   gene_id "CUFF.105"; transcript_id "CUFF.105.4"; FPKM "1.1485378983"; frac "0.233911"; conf_lo "0.853387"; conf_hi "1.443688"; cov "25.684547";
tig00000046 Cufflinks   transcript  29614   42766   181 +   .   gene_id "CUFF.105"; transcript_id "CUFF.105.5"; FPKM "0.2087076495"; frac "0.156792"; conf_lo "0.137112"; conf_hi "0.280304"; cov "4.667292";

What it came to my mind is:

sort the file by $4 (numeric: n) and $5 (reverse: nr) (sort -k4,4n -k5,5nr)
take the 1st line
remove unnecessary fields from $9

But this would fail for transcripts annotated in the reverse strand, and I won't be able to deal (how to fill) the $6  (score - a floating point value).
Any ideas of an existing method that could help me retrieve those gene lines?
Options to retrieve those lines directly from Cufflinks would also be accepted (I didn't find a suitable parameter in the manual).

Why is the score an issue? Can't you just keep the existing value? What am I missing?

Assuming none of the genes have duplicate entries across chromosomes (or strand or other biologically implausible things) in R one can:
library(rtracklayer)
library(GenomicRanges)
gtf = import.gff("file.gtf")
grl = split(gtf, gtf$gene_id)  # This produces a GRangesList
grl = endoapply(grl, function(x) {
    foo = x[1]
    foo$type = 'gene'
    start(foo) = min(start(x))
    end(foo) = max(end(x))
    return(c(foo, x))
    })
gr = unlist(gr)

You can then either write gr out to a file (either with export() or just printing the rows).

Typo on my part, I've updated that to import.gff().


text:
How to add OVERLAP_BP field to VEP output?
I'm using the Ensembl Variant Effect Predictor (VEP) to annotate VCF files, and I am trying to add some additional fields to the output, among these is OVERLAP_BP, which is listed here in available fields:

OVERLAP_BP - Number of base pairs overlapping with the corresponding variation feature

It's clear that some (many? all?) of these additional fields are only enabled by adding extra options to the VEP command (as listed on the RHS of the options table).
For example, the ZYG field is made available by adding the --individual flag to the VEP call.
But OVERLAP_BP doesn’t appear in the output fields column on that page, and when I add it to the output using --fields, it is consistently empty.
What options are required to make the OVERLAP_BP output field available to me in VEP?

I also opened this as a Github issue, and received the following response:

It seems there are some inconsistencies on our end. I'm looking into it.
  Until the fix, using --fields OverlapBP,OverlapPC works with any output format, while --fields OVERLAP_BP,OVERLAP_PC works only with default output format.
A small note that this overlap data is available for structural variants only.