text:
Protein sequence from patient data
Currently, I am working on NGS data and my aim is to get significance prediction of variants present in the vcf file. As we know about SIFT Score for significance score prediction, I am trying to understand how this score works. When I read its reference paper, in the method section, the procedure to determine the sift score starts with patient protein query sequence.
My question is how to get patient protein sequence from patient data (either vcf file or bam file) (which is not explained in the paper). Rest I can follow by other algorithms like PSI-BLAST & MOTIF.
Any suggestions....
Thanks

Ensembl now have a tool called Haplosaurus (still in beta) which will convert your phased VCF into actual protein sequences.

Haplosaurus
haplo is a local tool implementation of the same functionality that
  powers the Ensembl transcript haplotypes view. It takes phased
  genotypes from a VCF and constructs a pair of haplotype sequences for
  each overlapped transcript; these sequences are also translated into
  predicted protein haplotype sequences. Each variant haplotype sequence
  is aligned and compared to the reference, and an HGVS-like name is
  constructed representing its differences to the reference.
This approach offers an advantage over VEP's analysis, which treats
  each input variant independently. By considering the combined change
  contributed by all the variant alleles across a transcript, the
  compound effects the variants may have are correctly accounted for.
haplo shares much of the same command line functionality with vep, and
  can use VEP caches, Ensembl databases, GFF and GTF files as sources of
  transcript data; all vep command line flags relating to this
  functionality work the same with haplo.
Usage
Input data must be a VCF containing phased genotype data for at
  least one individual; no other formats are currently supported.
When using a VEP cache as the source of transcript annotation, the
  first time you run haplo with a particular cache it will spend some
  time scanning transcript locations in the cache.
./haplo -i input.vcf -o out.txt -cache

,@terdon : Thanks to both. I tried this. But I got some errors. Error is Attempt to reload Bio/EnsEMBL/VEP/TranscriptTree.pm aborted.. Although this file is present. I have checked  and confirmed this. I have also checked the similar issue at github (https://github.com/Ensembl/ensembl-vep/issues/169) which is still unsolved. Can you help me with this. Thanks.

You're best sticking with Irina on GitHub with this. She's working directly on developing the tool so will have the best insight and, if it's down to a bug, she can fix it.

I see that Irina has now solved this on GitHub. If anyone else has a similar problem, click through to the linked GitHub issue to see the resolution.

Thanks for software support. One last problem is stopping me to mark this as correct answer. I am not able to interpret output file. From your answer, I was expecting real protein sequence but I am getting many lines like this ENST00000342066 ENST00000342066:1027T>C ENSP00000342313:343W>R rs6672356 SM_5:2 in output file. Where is the protein sequence. How to read or interpret the output of ./haplo. Thanks.

The json option will give you sequence

Thanks for reply. I tried this too. The output file contains lots of sequences but it is very heavily packed with data or lots of data which is difficult to read. BTW, as you said, protein sequence can be obtained from phased vcf file which I got now. I have raised another issue related to performance of this tool please take a look at https://github.com/Ensembl/ensembl-vep/issues/214 . Thanks. Now I can mark it as right answer. Thanks again.

You can use SIFT in two different ways. One way is manually via their website. This is done by downloading the protein sequence (fasta) from a reliable repository (e.g., uniprot or NCBI), and change the amino acid of interest. Then calculate the SIFT score.
The other way, which is probably what you are interested in is with NGS data. Here normally SIFT is already known (precalculated) for all coding SNPs, and you can use that info directly during annotation of your variants. I don't know which variant annotation pipeline you are using, but for example annovar has SIFT data available at their databases repository.

You mean to say that for each variant at specific position (SNP especially), SIFT database has already determined which protein it will affect or what will be its corresponding protein sequence and how much similar it is with reference protein sequence and what will be its similarity score.  Is SIFT score valid for INDELs also.

Yes SIFT scores for the human genome are already precalculated, and you can use these scores for annotation of your variants (that is more efficient than to calculate them over and over again with each analysis, doesn't it?). SIFT works only for an amino acid substitution, that means nonsynonomous SNV. So no indels.


text:
Is there a standard definition for "assembly polishing"?

Is there a standard definition for "assembly polishing" in the field?
Is there a standard definition for what polishing algorithms do?

My understanding of "polishing" is strongly influenced by Pilon:
Whereby users input a reference FASTA along with the BAM of aligned reads against this reference FASTA. The polishing algorithm then does alignment to discover local misassemblies and other alignment inconsistencies. These are then "corrected", i.e. "polished".

Anything that makes something more correct (or shinier) would be polishing it in English, there's no special meaning here.

In long-read assembly, "polish" refers to the step to improve the base accuracy of contig sequences. I believe the terminology was originated from the HGAP paper:

The final consensus–calling algorithm Quiver, which takes into account all of the underlying data and the raw quality values inherent to SMRT sequencing, then polishes the assembly for final consensus accuracies in excess of 99.999% (QV of >50)...

Pilon was published more than a year later and in fact, neither its paper nor README mentions word "polish".
By the way, "error correction" refers to the step to fix sequencing errors in raw reads, not in contigs. In Canu/HGAP/Falcon, that step is done by pbdagcon or falcon_sense.

Disclaimer, this answer is based on feeling I got from talks and papers, but I do not have any hard reference supporting it.
I believe that genome polishing is a technique that was introduced for correction of individual bases of long noisy read assemblies - i.e. fixing SNPs and short indels. Guessing a correct base from long reads during assembly seems to be bit more complicated task in comparison to short read assemblies, therefore it is isolated into separated algorithms like Quiver/Arrow, Nanopolish (for correction by long reads) or the tool you mentioned - Pilon (correction by short reads).
Miss-assemblies could be corrected even prior long reads for instance by matching the assembly to optical mapping, or by tools like REAPR or TGnet. However, I have never heard term "genome polishing" used for this type of genome correction.

source: Stack Exchange

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

sort: votes


text:
Is there a way to do GWAS on phenotype data that is not normally distributed?
Is there a way to do GWAS on phenotype data that is not normally distributed?
For example, if you were measuring a trait as a proportion, and the majority of the data consisted of 0.00 with a long tail of variation. Or, perhaps the data was skewed the other way where the majority of the observations were 100% with a long tail towards zero.
Is there a GWAS method that is appropriate for these types of distributions?

Welcome to the site. I might be a bit dense but I don't understand the problem. In GWAS one looks which variants are associated with which trait. If you don't have information (or there is no variation in your dataset) for a given trait you can't relate the genome variants with it.

@Llopis There is still variation in a non-normal dataset. For example, 80% of your observations for a phenotypic trait may be measured at 0.00, but that leaves 20% at a value >0.00. Every test I have seen however requires the assumption that the data is normally distributed. My question is are there any methods for performing GWAS on non-normally distributed phenotype data.

Thanks for clarifying, now I understand better your problem. Perhaps it would help -to get an answer- to list some methods that assume normality in the data.

Yes, but like any other statistical test, you need to know the distribution of your data and pick an appropriate statistical test for it.

Yes, that can be done.
A common approach for these types of analysis is to carry out a transformation on the data first in order for it to have a normal distribution (the most general approach called an inverse normal transformation, among other names), then run the test on the transformed data.
As an example, for data that are positive and zero-skewed, a log transformation may be appropriate.
The transformation approach does require variation in the data values, though. For a trait that is exactly zero for most values, it might be better to treat it as a binary/categorical trait for the purpose of association testing, rather than a quantitative trait. Most association testing packages should allow associations for binary traits to be tested.

as mentioned in previous answers, transformations are frequently used. One commonly used method is quantile normal transformation. Basically you calculate the quantile from the original data and match it to a standard normal distribution. This transforms your data to a standard normal distribution. Even for traits that are naturally normal distributed like height, it is still a good idea to do this transformation as it has been shown to improve power as well.

Sure, you would use the negative binomial distribution when the majority of data was around 0. In other words the data is clustered with the standard deviation far greater than the mean. Most general linear modelling (GLM) will incorporate this approach.

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


text:
NCBI xml-format Printing nodeValues from tags with the same names
I would like to retrieve all Other-source_anchor on the level of DNA binding from within a Gene Ontology XML file using php.
If I try this with my code, I only get Gene Ontology back and not the nodeValues in Other-source_anchor. Could someone tell me how to do this?
Xml file were I wwant to get the nodeValues from
<Gene-commentary_heading>GeneOntology</Gene-commentary_heading>
      <Gene-commentary_source>
        <Other-source>
          <Other-source_pre-text>Provided by</Other-source_pre-text>
          <Other-source_anchor>GOA</Other-source_anchor>
          <Other-source_url>http://www.ebi.ac.uk/GOA/</Other-source_url>
        </Other-source>
      </Gene-commentary_source>
      <Gene-commentary_comment>
        <Gene-commentary>
          <Gene-commentary_type value="comment">254</Gene-commentary_type>
          <Gene-commentary_label>Function</Gene-commentary_label>
          <Gene-commentary_comment>
            <Gene-commentary>
              <Gene-commentary_type value="comment">254</Gene-commentary_type>
              <Gene-commentary_source>
                <Other-source>
                  <Other-source_src>
                    <Dbtag>
                      <Dbtag_db>GO</Dbtag_db>
                      <Dbtag_tag>
                        <Object-id>
                          <Object-id_id>3677</Object-id_id>
                        </Object-id>
                      </Dbtag_tag>
                    </Dbtag>
                  </Other-source_src>
                  <Other-source_anchor>DNA binding</Other-source_anchor>

code:
$esearch_test = "https://www.ncbi.nlm.nih.gov/gene/7128?report=xml&format=text";
$result = file_get_contents($esearch_test);
$xml = simplexml_load_string($result);
$doc = new DOMDocument();
$doc = DOMDocument::loadXML($xml);
$c = 1;
foreach($doc->getElementsByTagName('Gene-commentary_heading') as $node) {
    if ($node->textContent =="GeneOntology"){
        echo "
"."$c: ".$node->nodeName."
";
        echo "
"."$c: ".$node->nodeValue."
";
        $cn = 1;
        foreach ($node->getElementsByTagName('Other-source_anchor') as $subnode){
            echo "
"."$c: ".$subnode->nodeName."
";
            echo "
"."$c: ".$subnode->nodeValue."
";

            echo $cn++;
        }
    }
    $c++;
}

use a XSLT processor (http://php.net/manual/en/class.xsltprocessor.php) with the following stylesheet.
<?xml version='1.0' encoding="UTF-8"?>
<xsl:stylesheet xmlns:xsl='http://www.w3.org/1999/XSL/Transform' version='1.0'>
    <xsl:output method="html"/>

    <xsl:template match="/">

            <xsl:apply-templates select="//Gene-commentary_heading[ text() = 'GeneOntology']"/>

    </xsl:template>

    <xsl:template match="Gene-commentary_heading">

<xsl:value-of select="name(.)"/>

<xsl:value-of select="."/>

        <xsl:apply-templates select="..//Other-source_anchor"/>
    </xsl:template>

    <xsl:template match="Other-source_anchor">

<xsl:value-of select="name(.)"/>

<xsl:value-of select="."/>

    </xsl:template>

</xsl:stylesheet>


text:
What process and input data is required for a cellranger reference transcriptome?
I'm analysing single-cell RNA-Seq data using the 10X Genomics cellranger platform. While they provide reference data for Mouse and Humans, other species require a reference "transcriptome" to be generated with the cellranger mkref tool. However, I've found it to be very temperamental giving errors such as:
gffread genome_GFF3_genes.gff -T -o genome_genes.gtf

cellranger mkgtf genome_genes.gtf genome_genes.filtered.gtf \
                   --attribute=gene_biotype:protein_coding \
                   --attribute=gene_biotype:lincRNA \
                   --attribute=gene_biotype:antisense \
                   --attribute=gene_biotype:IG_LV_gene \
                   --attribute=gene_biotype:IG_V_gene \
                   --attribute=gene_biotype:IG_V_pseudogene \
                   --attribute=gene_biotype:IG_D_gene \
                   --attribute=gene_biotype:IG_J_gene \
                   --attribute=gene_biotype:IG_J_pseudogene \
                   --attribute=gene_biotype:IG_C_gene \
                   --attribute=gene_biotype:IG_C_pseudogene \
                   --attribute=gene_biotype:TR_V_gene \
                   --attribute=gene_biotype:TR_V_pseudogene \
                   --attribute=gene_biotype:TR_D_gene \
                   --attribute=gene_biotype:TR_J_gene \
                   --attribute=gene_biotype:TR_J_pseudogene \
                   --attribute=gene_biotype:TR_C_gene

cellranger mkref --genome="reference-name" \
                 --fasta="transcriptome.fasta" \
                 --genes="genome_genes.filtered.gtf"

Fatal INPUT FILE error, no valid exon lines in the GTF file: reference_name/genes/genes.gtf
  Solution: check the formatting of the GTF file. Most likely cause is the difference in chromosome naming between GTF and FASTA file.
subprocess.CalledProcessError: Command '['STAR', '--runMode', 'genomeGenerate', '--genomeDir', 'reference_name/star', '--runThreadN', '1', '--genomeFastaFiles', 'reference_name/fasta/genome.fa', '--sjdbGTFfile', 'reference_name/genes/genes.gtf', '--limitGenomeGenerateRAM', '17179869184', '--genomeSAsparseD', '1', '--genomeSAindexNbases', '12', '--genomeChrBinNbits', '10']' returned non-zero exit status 104
subprocess.CalledProcessError: Command '['STAR', '--runMode', 'genomeGenerate', '--genomeDir', 'reference_name/star', '--runThreadN', '1', '--genomeFastaFiles', 'reference_name/fasta/genome.fa', '--sjdbGTFfile', 'reference_name/genes/genes.gtf', '--limitGenomeGenerateRAM', '17179869184', '--genomeSAsparseD', '1', '--genomeSAindexNbases', '12', '--genomeChrBinNbits', '10']' returned non-zero exit status -11

What input data is required for the cellranger to construct a reference? What format do the fasta and gtf files need be in? Is it possible to use a "transcriptome" rather than a "genome" as a reference here? For example, I have a transcriptome reference with fasta names >name for each transcript rather than each chromosome.

Your problem is caused by using the transcriptome fasta file rather than the genome fasta file. You've already given it transcriptome information with genome_genes.filtered.gtf, it needs the genomic sequence to go along with that.
As an aside, your GTF filtering command is unlikely to be useful unless you have a Gencode or Ensembl GTF for human. Have a look at the gene_biotypes in your file and decide for yourself what (if anything) needs to be filtered. For example, I made a cellranger reference last night for dm3 (the old fruit fly genome) with:
cellranger mkgtf annotation.gtf genes.gtf \
    --attribute=gene_biotype:protein_coding \
    --attribute=gene_biotype:ncRNA \
    --attribute=gene_biotype:snRNA \
    --attribute=gene_biotype:snoRNA \
    --attribute=gene_biotype:pre_miRNA

That should be customized for relevance. It's quite possible that your GTF doesn't even have gene_biotype attributes, that's specific to Ensembl I think.

Indeed, I am not working with Humans on this project (I tried to include a "model-organisms" tag to emphasise this). I included this from the cellranger documentation while troubleshooting the issue. To clarify, the reference "genome" is a FASTA with chromosomes named? I have a recent annotation (GTF) for A. thaliana but the reference (FASTA) is by transcripts in this database. An older full-genome reference is available but I'd rather use a recently updated on is possible.

Yes, the reference genome is the fasta file containing chromosomal entries. You can get appropriate GTF files from Ensembl, a transcriptome fasta file is of no use to you.


text:
verifybamid results with freemix < 2% and chipmix > 10%
I ran verifyBAMID on each 200 chip-seq BAM files with input merged vcf file of Genotypes of 600 samples ( 200 of them are chip-seq individuals). I had AF,AC,AN in the info column which was calculated using these genotypes.
~80% of 200 samples are Caucasians.
I am struggling to decide for the samples which have freemix < 2% but chipmix >2%.
For sample 1365, I had histone modified chip-seq sample 1365_A with two input samples 1365_Pos_IN_2, 1365_Pos_IN_1.

histone modified chip-seq BAM 1365_A has .00669  whereas chipmix is 0.3765 where with input samples it does resonable job with correct matching of genotype column and reasonable chipmix and freemix values.
Why is chipmix so high for histone sample?
SEQ_ID  RG  CHIP_ID #SNPS   #READS  AVG_DP  FREEMIX FREELK1 FREELK0 FREE_RH FREE_RA CHIPMIX CHIPLK1 CHIPLK0 CHIP_RH CHIP_RA DPREF   RDPHET  RDPALT
1365_A  ALL 1365_A  683278  436495  0.64    **0.00669** 160631.69   160867.97   NA  NA  **0.37645** 145839.46   194576.92   NA  NA  0.657   0.9478  0.9369
1365_Pos_IN_2   ALL 1365_A  683278  1770725 2.59    **0.00320** 676736.66   676773.49   NA  NA  **0.00298** 425582.72   426459.54   NA  NA  2.598   0.9965  0.9923
1365_Pos_IN_1   ALL 1365_A  683278  1127574 1.65    **0.00000** 446421.43   446421.43   NA  NA  **0.00193** 254020.46   254353.15   NA  NA  1.653   0.9961  0.9982

What does it mean when your target sample matched with the intended sample genotype from a vcf file with < 2% free mix value, but with chipmix value > 2%?
for sample 1368_A and 1213_A
1368_A  ALL 1368_A  683278  450800  0.66    **0.01019** 167252.18   167673.56   NA  NA  **0.56004** 158591.79   248605.77   NA  NA  0.677   0.9507  0.9477
1213_A  ALL 1213_A  683278  945562  1.38    **0.01361** 345881.30   347208.39   NA  NA  **0.97533** 346470.85   647079.17   NA  NA  1.396   0.9783  0.9891

This is clearly contaminated sample with 1287_A. then why is its chipmix not 0.9?
1215_Pos_IN ALL 1287_A  683278  3336277 4.88    **0.15412** 1360018.71  1407231.58  NA  NA  **0.18590** 1109098.83  1431053.16  NA  NA  4.900   0.9941  0.9917


text:
Automating download of results of Drug and Drug Target Mapping
The drug and drug target mapping tool shows the top 3 PDB IDs for each drug based on the drug target sequence similarity search and ligand id search. Is there a way to get all the PDB ids of the drug-target PDB files and not only the top 3?
I know you can get this information after clicking the "find" button under the last column of the output table (the column named "All PDB"), but I have a query of several hundreds of drugs and I'm not willing to do it manually for each result.

Simpler option would be to not use the query interface  and directly parse the CSV file found on that page.
I just read @llrs comment and I guess it's very much the same spirit as my own - download the data and query it programmatically. I see now (was on mobile earlier) that the 'FIND' URL links out to a BLAST search. so it might be as easy as retrieving the relevant protein sequence and replicating the URL. Or if you want to go a different route and finding all PDB identifiers with the given ligand, cross-reference the above data with what you'd find here


text:
bcftools consensus not working when using -H flag
Sorry, I have been trying to work this out desparately all day. Biostars couldn't find any answer so I am cross posting here.
I am trying to get both the left and right reference genome using a sample name in phase VCF file.
This works:
bcftools consensus -f lyrata_chr01-short.fasta phasedVCF-short.vcf.gz -s ms01e -H 1 > ms01e-left.fa

When I am trying to make the new Ref Sequence by imputing genotype from the right haplotype, this doesn't work.
bcftools consensus -f lyrata_chr01-short.fasta phasedVCF-short.vcf.gz -s ms01e -H 2 > ms01e-right.fa
Broken VCF, too few alts at 1:141

What is the issue ?
Here is my VCF data:
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  MA605   MA611   MA622   MA625   MA629   Ncm8    Sp154   Sp164   Sp21    Sp3 Sp76    SpNor33 ms01e   ms02g   ms03g   ms04h
1   141 .   C   T   .   .   .   GT:PG:PG_al:PI  .   .   .   0|1:0|1:C|T:10763   .   .   .   .   .   .   .   .   .   .   .   .
1   3697    .   T   C   .   .   .   GT:PG:PG_al:PI  .   0|0:0|0:T|T:11398   .   0|1:0|1:T|C:8039    0|1:0|1:T|C:8202    0|1:0|1:T|C:4616    0|1:0|1:T|C:3938    .   0|1:0|1:T|C:5102    .   .   0|1:0|1:T|C:17340   1|0:1|0:C|T:1   .   .   .

The GT field is empty at POS 141 for SAMPLE ms01e, but in that situation it should by default pick a REF allele. Isn't it?
And, it only is the problem when doing -H 2 not -H 1.
Update: I have been trying to fix this issue by playing with -i and -e using the BCFtools manual. But, still cannot find the fix.

I don't know much about VCF but the first row seems to have less fields than the second. Is this normal/accepted format ? Also did you try the latest version (as per ) ?


text:
Viral Load and Antiretroviral Therapy Dataset
My project is on Mathematical Modeling of the relationship between viral load and antiretroviral therapy. Its purely an academic investigation which is exploratory in nature.
Question I am seeking a public-database or a published data set for drug-HIV interactions and in particular drug vs viral load, for example a time series.
Project Intend using the different orthogonal polynomials to fit the curve and then test them statistically to know which fits best.
Any other suggestions are welcomed.

Which curve do you want to fit? To which data? Or are you asking us to where to find data of antiretrovirals and viral load?


text:
how to do analysis of a table content SNP?
I have a table with information about SNPs, the table has a total of 55127 rows.
How can I know the number of gene mutations in the table and if a gene is present?
The table is in xls format.

Hi Sofia, recommend  a package for what? Could you [edit] the question and explain what is your goal? Also I assume you mean SNP or what is a SNIP? As a side note all post are already signed with the account of the creator, please do not sign them in the body of the post.

Thanks for the edits, but the main question stills remain. What is your goal, what do you want to do with this data? Do you have a question in mind, a disease, a problem, a treatment? Could you be more clear on this regard, otherwise is too unclear to try an answer

thanks for the answer, i want to know the number of gene mutations. and to know if a certain gene is in the table.

It would be helpful to know the format of the data and the file size. Can you post an example of what the data look like.

Please [edit] your question and explain exactly what you are trying to do. What output do you need? Do you want a list of SNPs per gene? A list of genes in the file? None of these need to be done in R (nor is R particularly well suited for this sort of simple text parsing), do you need R? Also, what operating system are you using? We need to know that to know what tools you have easily available. Finally, what format is the file in? Is that an excel file? A csv? Tab separated? Space separated?

For working with a relatively small dataset in tabular format, I would recommend dplyr. You can read in a excel file using readxl
library(readxl)
library(dplyr)

df <- read_excel(FILENAME)

To get the number of mutations per gene:
df %>% group_by(`Gene name`) %>%
  summarise(num_mutations=n())

If you want to know if there are mutations in a particular gene:
df %>% filter(`Gene name` == QUERY)


text:
Transform a R data.frame into a krona plot without krona tools?
KronaTools lacks a R pipeline to convert any data.frame into a krona plot. There is a way to go via a phyloseq-class object, but coming from a data.frame I did not find any documentation on how to transform these.

First download any Krona plot (.html) file or make one by using the Excel sheet from the Krona Package. Open in text-editor and remove anything below "datasets" tag and/or the first "node" tag till the end of the document.
Here is the path to go: 1) get your data.frame in hierachical order 2) transform to data.tree object 3) transform to xml file 4) gsub until all the tags match krona annotation
options(scipen=999)
options(stringsAsFactors=FALSE)
library(data.tree)
library(xml2)
library(magrittr)

#get your data into a hieracial format in a data.frame
t <- structure(list(freq = c(34L, 127L, 42L, 123L, 343L, 13L, 104L,
21L, 135L, 91L, 167L, 148L, 90L, 46L, 25L, 79L, 1726L), king.n = c("",
"Fungi", "Fungi", "Fungi", "Fungi", "Fungi", "Fungi", "Plantae",
"Plantae", "Plantae", "Plantae", "Plantae", "Plantae", "Plantae",
"Plantae", "Plantae", "Plantae"), class.n = c("", "Dothideomycetes",
"Dothideomycetes", "Dothideomycetes", "Dothideomycetes",     "Dothideomycetes",
"Dothideomycetes", "Liliopsida", "Magnoliopsida", "Magnoliopsida",
"Magnoliopsida", "Magnoliopsida", "Magnoliopsida", "Magnoliopsida",
"Magnoliopsida", "Magnoliopsida", "Magnoliopsida"), fam.n = c("",
"", "Cladosporiaceae", "Didymellaceae", "Pleosporaceae", "Pleosporaceae",
"Pleosporaceae", "Poaceae", "", "Asteraceae", "Asteraceae", "Asteraceae",
"Asteraceae", "Fabaceae", "Hydrophyllaceae", "Hypericaceae",
"Malvaceae"), gen.n = c("", "", "Cladosporium", "", "Alternaria",
"Pyrenophora", "Pyrenophora", "Phleum", "", "", "Achillea", "Centaurea",
"Jurinea", "Trifolium", "Phacelia", "Hypericum", "Malva"), spec.n = c("",
"", "cladosporioides", "", "", "", "tritici-repentis", "pratense",
"", "", "", "cyanus", "", "", "", "", "")), .Names = c("freq",
"king.n", "class.n", "fam.n", "gen.n", "spec.n"), row.names = c(1L,
3L, 4L, 5L, 8L, 9L, 11L, 15L, 16L, 17L, 18L, 19L, 21L, 22L, 25L,
26L, 28L), class = "data.frame")

#if any sequences can not be assigned to the highest hieracial level give them a random name like "NA"
t[t[,"king.n"]=="","king.n"] <- "NA"
#define the hieracy (do not include the counting var here), the first space is the root category
t$pathString <- paste(" ",t$king.n,t$class.n,t$fam.n,t$gen.n,t$spec.n,sep = "/")
#transform to a data.tree
t <- as.Node(t)
#define the sequence frequency calculations
t$Do(function(x) x$freq <- Aggregate(x, "freq", sum), traversal = "post-order")
print(t,"freq")
#get all frequencies for each node
repval <- t$Get("freq")

#tranform into a nested list and afterwards into a xml document to get the nesting
t <- t %>% as.list(.,mode="simple",unname=T)
t <- as_xml_document(list(root=t))%>%as.character

#from here it is just a bunch of regex
#split the xml file into one line for each <
temp <- strsplit(t,"\n")[[1]][2] %>% strsplit(.,"(?=<)",perl=T) %>% .[[1]]
temp <- paste0(temp[c(T,F)],temp[c(F,T)]) %>% gsub("[0-9]","",.)
#remove the names nodes and replace by krona annotation
#replace the closing named tags by generic closing "node" tags
temp <- gsub("</(.*)?>","</node>",temp) %>% gsub(">.*$",">",.)
#replace the named opening tags by generic "node" tags with the "name" attribute, replace all (because they are incomplete) numeric values with a placeholder
temp <- gsub("<(?=[^/])(.*)?>",
     '<node name="\\1"><magnitude><val>__\\1</val></magnitude><score><val>1</val></score>',temp,perl=T)
temp %<>% gsub("root"," ",.)
#find the placeholders to insert the frequencies and execute
repos <- grep("__",temp)
for (i in repos)
  temp[i] <- gsub("^(.*)(__.*?)(?=<)(.*)",paste0("\\1",repval[which(repos==i)],"\\3"),temp[i],perl=T)

#add the end of the .html file
temp <- c(temp,"</krona>
","","")  %>% unlist
tout <- paste(temp,collapse="\n")

#append your data to the Krona plot
write(tout,file="Krona1.html",append=TRUE)

Et voilá!

Did you just question and answer your own question at the same time?

@Shique that is even encouraged: Can I answer my own question?

@Jaap I am aware of that. However I thought it was meant to be used if no one could answer it in a certain amount of time. More like "While I was trying to fix this on my own, I found out that... " Rather than "Hey I just found something, let me post this on SO as a how-to"


text:
How to compare the contents of a column of the same data table
I have this table:

and I want to get the rows that are equal to the first three columns, like this:

I've tried these functions, but when I get the index of the lines, r doesn't give the output that I want:
df$obj<-sapply(c("sample1", "sample2", "sample3"), function (comparison) {
    df$sample1 == df[comparison]
})

following on from my dplyr suggestion on your previous question: filter(df, sample1 == sample2 & sample1 == sample3) Note that a and A are not the same. If you want a case-insensitive matching, use filter(df, toupper(sample1) == toupper(sample2) & toupper(sample1) == toupper(sample3))

but in that way you only get comparison between two columns, you see? you can have the same result between sample 1 and sample 2 , but when you are going to compare sample 1 and sample 3, the result can be different than the previous one

I'm not sure what you mean. If sample1 == sample2 and sample1 == sample3, then sample2 has to be equal to sample3. Unless I'm misunderstanding your goal. My code will turn your example input into your desired output though

Loop through 3 columns by row using apply, then convert toupper case (as we want a same as A), get unique, and get lengths, if length is 1 then all 3 values are equal:
df1[ lengths(
  apply(df1[, 1:3], 1, function(i) unique(toupper(i)))
  ) == 1, ]

This would scale better if we want to compare 4, 10, n(?) columns.

No need to use dplyr, we can do this using base, see below example:
# reproducible example data
df1 <- data.frame(
  t(data.frame(
    one = c("A", "A", "A", "C", "C", "C"),
    two = c("A", "A", "A", "A", "A", "C"),
    thr = c("A", "A", "A", "A", "A", "A"),
    fou = c("C", "A", "A", "A", "A", "A"))),
  stringsAsFactors = FALSE)

df1
#     X1 X2 X3 X4 X5 X6
# one  A  A  A  C  C  C
# two  A  A  A  A  A  C
# thr  A  A  A  A  A  A
# fou  C  A  A  A  A  A

# solution, subset using logical index:
selection <- df1[df1[, 1] == df1[, 2] & df1[, 1] == df1[, 3], ]

selection
#     X1 X2 X3 X4 X5 X6
# one  A  A  A  C  C  C
# two  A  A  A  A  A  C
# thr  A  A  A  A  A  A

To be fair, you don't really need dplyr for anything. I just think it's easier to wrap your head around filter then a command that repeats the name of the data frame 5 times, along with 5 sets of square brackets and 5 commas. Personal preference though I guess.

@heathobrien, you are right dplyr is a great tool, and I think for some difficult tasks even necessary. For row selection with conditional statements I am used to select this (old fashion) way. But with dplyr would certainly be good as well.

t(data.frame(…)) is an odd pattern, since t always converts data.frames to matrices, and assigning the result to a variable called df is actively misleading.

@KonradRudolph, that was to make a workable example. I understand that you don't like my style (also from previous comments), feel free to put your own answer.

@b.nota Actually I have no problem with the style of this answer, I think it actually works well otherwise. But, to clarify, my previous comment wasn’t about style.

@b.nota hope you don't mind the edit. Generally I avoid function names (df). And converted your matrix to a dataframe by wrapping it with another dataframe. Your solution works just fine.


text:
Why isn't there a standard method to convert GFA to JSON?
I'm confused why the GFA format doesn't have a standard JSON parser. It appears to be related to a decision made several years ago, here.
The objection appears to be "JSON isn't necessary, as it's just a simple linked list. Users could write their own". However, I would think this is a useful feature.
There are applications I could think of whereby a JSON would be necessary, or at least very convenient.
(1) Is there a standard way to convert GFA into JSON? Or do users write their own scripts?
(2) Is there a conceptual reason why GFA doesn't convert into JSON?

I thought the explanation Heng gave was pretty comprehensive: it’s overkill, and you can use JSON in tags. Of course you shouldn’t write your own parser, there are enough existing ones.

@KonradRudolph "there are enough existing ones." Do you have recommendations? Heng Li mentions his own in that thread, but there are others?

My recommendation is to use the standard library in whatever language you’re using. If that language doesn’t have a JSON parser in the standard library, use the most widely used implementation.

@KonradRudolph "use the most widely used implementation" What is this?

That’s the point: it depends on the programming language. A simple search of “JSON” + name of the programming language should have it as the first or second hit.

@KonradRudolph Oh, I understand what you meant by the comment above now. I was thinking you meant a specific JSON parser for GFAs.

Short answer: no, there isn't a "standard" method, and there won't be.
Long answer: like TAB-delimited and XML, JSON alone is not a specific format. You have to define a schema to give meanings to data. However, unlike TAB-delimited "format" and XML, there is not an official way to define JSON schema. There have been various attempts but all of them are complicated and none of them are widely adopted. As you can't write down a JSON schema for GFA in the first place, there won't be a "standard" way to convert GFA to JSON. The best you can do is to roll out your own in-house solution, unfortunately.

Thanks. "There have been various attempts but all of them are complicated and none of them are widely adopted." Just to clarify though, why couldn't the community adopt a JSON scheme for GFA? If you encode all information of a GFA into one JSON format, it would be easy for others to modify this if they wish, right?

You are now asking a different question: why not use JSON to encode graphs? This question was answered in the thread you quoted.

" in memory, a JSON object is also fatter, which can be a concern for huge graphs." This reason?


text:
GFF file too big to load into R
I want to load a GFF (annotation) file into an R data frame in order to extract some information from it, e.g. the location of transposons in the genome.
I used the following code:
library(rtracklayer)
gff <- rtracklayer::import('c_elegans.PRJNA13758.WS263.annotations.gff3.gz')
Error: cannot allocate vector of size 210.2 Mb

It seems that the GFF3 file is too big for loading into R. Do you know if this is limited by R or by my computer's memory?
I ran this on a GNU/Linux computer.

It depends of course on your available RAM on your computer.
Under Windows you can increase the available RAM in R with memory limit:
> memory.limit()
[1] 8070
> memory.limit(size=10000)
[1] 10000
> memory.limit()
[1] 10000

But that is not necessary on linux/unix.

This only works (and is only required) on Windows. OP is using Linux. Not sure why he accepted this answer.

Indeed, memory.limit()
[1] Inf
Warning message:
'memory.limit()' is Windows-specific. I understood this as 'no memory limit' on Linux, so only limited by available RAM?


text:
single stranded cDNA vs double-stranded cDNA
Single stranded cDNA is used by race pcr and rt-pcr while next-gen sequencing (Illumina) requires double stranded cDNA. What is the reason behind this ?

I'm voting to close this question as off-topic because the question is not about bioinformatics, but for some specific methods of the wet lab. If you have a question about how to handle these sequences by computers, I think it would be better suited for the site.

Can we migrate this to biology? It's a valid question and I think it's even somewhat relevant to bioinformatic (it's important to understand generation of the data), but it would suit better to Biology.SE I think. Anyway, could you link source of your information?

Sequencing also requires single stranded DNA, otherwise the sequencing primers can’t bind. Sure there are steps in the protocol where double stranded material exists, but that’s transient.


text:
Plotting gene models - alternatives to ggbio
I am curious how people make genome track plots such as what can be done with the Bioconductor package ggbio.
For example this plot was made w/ ggbio:

I think ggbio is great for making fast plots, but it is kind of a pain to make publication quality figures because it is hard to manipulate the plots after they are constructed. I have noticed that ggbio constructed plots do not return a ggplot or gtable object. I can access the grobs from a ggbio tracks object, but they won't be aligned as they would be in the final plot. I like to make shared legends and include them in the plot, and without having a gtable or ggplot object it is hard to add the shared legends in, and I end up having to save them as a separate file which is tedious. FYI I mainly use ggbio for gene tracks + x-axis alignment.
Normally I would use ggplot2 to customize the plot extensively, but it doesn't work very well after plotting with ggbio.
How can I use ggbio for these types of plots? Is there a better tool instead?
EDIT:
I would get the legends for each individual plot and combine them before adding back to the main plot using something like this:
leg1.grob = cowplot::get_legend(plot1)
leg2.grob = cowplot::get_legend(plot2)
sharedleg.grob = gridExtra::gtable_rbind(leg1,leg2)

mainplot.grob = ggplotGrob(main_plot)
mainplot.grob = gtable_add_cols(mainplot.grob,
                                widths = sum(sharedleg.grob$widths))
mainplot.grob = gtable_add_grobs(mainplot.grob,
                                 grobs = sharedleg.grob,
                                 t=t,b=b,l=l,r=r)
cowplot::ggdraw(mainplot.grob)

I've used gviz before, but long ago. Do you know if you can extract the grid plot grobs of the gene track plots with Gviz?

It seems like the ggbio object already contains the corresponding ggplot object.
You should be able to extract the ggplot from the ggbio_obj@ggplot slot and extract your gtables.
Small example:
library(ggbio)
library(ggplot2)
library(patchwork)

fl.bam <- system.file("extdata", "wg-brca1.sorted.bam", package = "biovizBase")
wh <- as(c("chr17:41239394-41319151:+"), "GRanges")
p_cov <- autoplot(fl.bam, which = wh)
#> reading in as Bamfile
#> Parsing raw coverage...
#> Read GAlignments from BamFile...
#> extracting information...

df <-  data.frame(x= sample(1:100,size = 100,replace = TRUE))
p_hist <- ggplot(df) +
    geom_histogram(aes(x))

p_hist / p_cov
#> Error: Can't add `e2` to a ggplot object.
#> Backtrace:
#>     █
#>  1. └─patchwork:::`/.ggplot`(p_hist, p_cov)
#>  2.   └─GGally:::`+.gg`(e1, e2)
#>  3.     └─e1 %+% e2
#>  4.       └─ggplot2:::add_ggplot(e1, e2, e2name)
#>  5.         ├─ggplot2::ggplot_add(object, p, objectname)
#>  6.         └─ggplot2:::ggplot_add.default(object, p, objectname)

p_hist / p_cov@ggplot
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Created on 2021-02-05 by the reprex package (v1.0.0)
Session info

sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R version 3.6.0 (2019-04-26)
#>  os       Red Hat Enterprise Linux
#>  system   x86_64, linux-gnu
#>  ui       X11
#>  language (EN)
#>  collate  en_US.UTF-8
#>  ctype    en_US.UTF-8
#>  tz       Europe/Berlin
#>  date     2021-02-05
#>
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package              * version  date       lib source
#>  acepack                1.4.1    2016-10-29 [1] CRAN (R 3.6.0)
#>  AnnotationDbi          1.48.0   2019-10-29 [1] Bioconductor
#>  AnnotationFilter       1.10.0   2019-10-29 [1] Bioconductor
#>  askpass                1.1      2019-01-13 [1] CRAN (R 3.6.0)
#>  assertthat             0.2.1    2019-03-21 [1] CRAN (R 3.6.0)
#>  backports              1.1.8    2020-06-17 [1] CRAN (R 3.6.0)
#>  base64enc              0.1-3    2015-07-28 [1] CRAN (R 3.6.0)
#>  Biobase                2.46.0   2019-10-29 [1] Bioconductor
#>  BiocFileCache          1.10.2   2019-11-08 [1] Bioconductor
#>  BiocGenerics         * 0.32.0   2019-10-29 [1] Bioconductor
#>  BiocManager            1.30.10  2019-11-16 [1] CRAN (R 3.6.0)
#>  BiocParallel           1.20.1   2019-12-21 [1] Bioconductor
#>  biomaRt                2.42.1   2020-03-26 [1] Bioconductor
#>  Biostrings             2.54.0   2019-10-29 [1] Bioconductor
#>  biovizBase             1.34.1   2019-12-04 [1] Bioconductor
#>  bit                    1.1-15.2 2020-02-10 [1] CRAN (R 3.6.0)
#>  bit64                  0.9-7    2017-05-08 [1] CRAN (R 3.6.0)
#>  bitops                 1.0-6    2013-08-17 [1] CRAN (R 3.6.0)
#>  blob                   1.2.1    2020-01-20 [1] CRAN (R 3.6.0)
#>  BSgenome               1.54.0   2019-10-29 [1] Bioconductor
#>  checkmate              2.0.0    2020-02-06 [1] CRAN (R 3.6.0)
#>  cli                    2.0.2    2020-02-28 [1] CRAN (R 3.6.0)
#>  cluster                2.0.8    2019-04-05 [2] CRAN (R 3.6.0)
#>  colorspace             1.4-1    2019-03-18 [1] CRAN (R 3.6.0)
#>  crayon                 1.3.4    2017-09-16 [1] CRAN (R 3.6.0)
#>  curl                   4.3      2019-12-02 [1] CRAN (R 3.6.0)
#>  data.table             1.12.8   2019-12-09 [1] CRAN (R 3.6.0)
#>  DBI                    1.1.0    2019-12-15 [1] CRAN (R 3.6.0)
#>  dbplyr                 1.4.2    2019-06-17 [1] CRAN (R 3.6.0)
#>  DelayedArray           0.12.2   2020-01-06 [1] Bioconductor
#>  dichromat              2.0-0    2013-01-24 [1] CRAN (R 3.6.0)
#>  digest                 0.6.25   2020-02-23 [1] CRAN (R 3.6.0)
#>  dplyr                  0.8.5    2020-03-07 [1] CRAN (R 3.6.0)
#>  ellipsis               0.3.1    2020-05-15 [1] CRAN (R 3.6.0)
#>  ensembldb              2.10.2   2019-11-20 [1] Bioconductor
#>  evaluate               0.14     2019-05-28 [1] CRAN (R 3.6.0)
#>  fansi                  0.4.1    2020-01-08 [1] CRAN (R 3.6.0)
#>  farver                 2.0.3    2020-01-16 [1] CRAN (R 3.6.0)
#>  foreign                0.8-71   2018-07-20 [2] CRAN (R 3.6.0)
#>  Formula                1.2-3    2018-05-03 [1] CRAN (R 3.6.0)
#>  fs                     1.4.0    2020-03-31 [1] CRAN (R 3.6.0)
#>  GenomeInfoDb           1.22.1   2020-03-27 [1] Bioconductor
#>  GenomeInfoDbData       1.2.2    2020-04-01 [1] Bioconductor
#>  GenomicAlignments      1.22.1   2019-11-12 [1] Bioconductor
#>  GenomicFeatures        1.38.2   2020-02-15 [1] Bioconductor
#>  GenomicRanges          1.38.0   2019-10-29 [1] Bioconductor
#>  GGally                 1.5.0    2020-03-25 [1] CRAN (R 3.6.0)
#>  ggbio                * 1.34.0   2019-10-29 [1] Bioconductor
#>  ggplot2              * 3.3.2    2020-06-19 [1] CRAN (R 3.6.0)
#>  glue                   1.4.1    2020-05-13 [1] CRAN (R 3.6.0)
#>  graph                  1.64.0   2019-10-29 [1] Bioconductor
#>  gridExtra              2.3      2017-09-09 [1] CRAN (R 3.6.0)
#>  gtable                 0.3.0    2019-03-25 [1] CRAN (R 3.6.0)
#>  highr                  0.8      2019-03-20 [1] CRAN (R 3.6.0)
#>  Hmisc                  4.4-0    2020-03-23 [1] CRAN (R 3.6.0)
#>  hms                    0.5.3    2020-01-08 [1] CRAN (R 3.6.0)
#>  htmlTable              1.13.3   2019-12-04 [1] CRAN (R 3.6.0)
#>  htmltools              0.5.0    2020-06-16 [1] CRAN (R 3.6.0)
#>  htmlwidgets            1.5.1    2019-10-08 [1] CRAN (R 3.6.0)
#>  httr                   1.4.1    2019-08-05 [1] CRAN (R 3.6.0)
#>  IRanges                2.20.2   2020-01-13 [1] Bioconductor
#>  jpeg                   0.1-8.1  2019-10-24 [1] CRAN (R 3.6.0)
#>  knitr                  1.29     2020-06-23 [1] CRAN (R 3.6.0)
#>  labeling               0.3      2014-08-23 [1] CRAN (R 3.6.0)
#>  lattice                0.20-38  2018-11-04 [2] CRAN (R 3.6.0)
#>  latticeExtra           0.6-29   2019-12-19 [1] CRAN (R 3.6.0)
#>  lazyeval               0.2.2    2019-03-15 [1] CRAN (R 3.6.0)
#>  lifecycle              0.2.0    2020-03-06 [1] CRAN (R 3.6.0)
#>  magrittr               1.5      2014-11-22 [1] CRAN (R 3.6.0)
#>  Matrix                 1.2-17   2019-03-22 [2] CRAN (R 3.6.0)
#>  matrixStats            0.56.0   2020-03-13 [1] CRAN (R 3.6.0)
#>  memoise                1.1.0    2017-04-21 [1] CRAN (R 3.6.0)
#>  mime                   0.9      2020-02-04 [1] CRAN (R 3.6.0)
#>  munsell                0.5.0    2018-06-12 [1] CRAN (R 3.6.0)
#>  nnet                   7.3-12   2016-02-02 [2] CRAN (R 3.6.0)
#>  openssl                1.4.1    2019-07-18 [1] CRAN (R 3.6.0)
#>  OrganismDbi            1.28.0   2019-10-29 [1] Bioconductor
#>  patchwork            * 1.0.0    2019-12-01 [1] CRAN (R 3.6.0)
#>  pillar                 1.4.3    2019-12-20 [1] CRAN (R 3.6.0)
#>  pkgconfig              2.0.3    2019-09-22 [1] CRAN (R 3.6.0)
#>  plyr                   1.8.6    2020-03-03 [1] CRAN (R 3.6.0)
#>  png                    0.1-7    2013-12-03 [1] CRAN (R 3.6.0)
#>  prettyunits            1.1.1    2020-01-24 [1] CRAN (R 3.6.0)
#>  progress               1.2.2    2019-05-16 [1] CRAN (R 3.6.0)
#>  ProtGenerics           1.18.0   2019-10-29 [1] Bioconductor
#>  purrr                  0.3.4    2020-04-17 [1] CRAN (R 3.6.0)
#>  R6                     2.4.1    2019-11-12 [1] CRAN (R 3.6.0)
#>  rappdirs               0.3.1    2016-03-28 [1] CRAN (R 3.6.0)
#>  RBGL                   1.62.1   2019-10-30 [1] Bioconductor
#>  RColorBrewer           1.1-2    2014-12-07 [1] CRAN (R 3.6.0)
#>  Rcpp                   1.0.4.6  2020-04-09 [1] CRAN (R 3.6.0)
#>  RCurl                  1.98-1.1 2020-01-19 [1] CRAN (R 3.6.0)
#>  reprex                 1.0.0    2021-01-27 [1] CRAN (R 3.6.0)
#>  reshape                0.8.8    2018-10-23 [1] CRAN (R 3.6.0)
#>  reshape2               1.4.4    2020-04-09 [1] CRAN (R 3.6.0)
#>  rlang                  0.4.6    2020-05-02 [1] CRAN (R 3.6.0)
#>  rmarkdown              2.1      2020-01-20 [1] CRAN (R 3.6.0)
#>  rpart                  4.1-15   2019-04-12 [2] CRAN (R 3.6.0)
#>  Rsamtools              2.2.3    2020-02-23 [1] Bioconductor
#>  RSQLite                2.2.0    2020-01-07 [1] CRAN (R 3.6.0)
#>  rstudioapi             0.11     2020-02-07 [1] CRAN (R 3.6.0)
#>  rtracklayer            1.46.0   2019-10-29 [1] Bioconductor
#>  S4Vectors              0.24.3   2020-01-18 [1] Bioconductor
#>  scales                 1.1.1    2020-05-11 [1] CRAN (R 3.6.0)
#>  sessioninfo            1.1.1    2018-11-05 [1] CRAN (R 3.6.0)
#>  stringi                1.4.6    2020-02-17 [1] CRAN (R 3.6.0)
#>  stringr                1.4.0    2019-02-10 [1] CRAN (R 3.6.0)
#>  styler                 1.3.2    2020-02-23 [1] CRAN (R 3.6.0)
#>  SummarizedExperiment   1.16.1   2019-12-19 [1] Bioconductor
#>  survival               3.1-11   2020-03-07 [1] CRAN (R 3.6.0)
#>  tibble                 3.0.0    2020-03-30 [1] CRAN (R 3.6.0)
#>  tidyselect             1.0.0    2020-01-27 [1] CRAN (R 3.6.0)
#>  VariantAnnotation      1.32.0   2019-10-29 [1] Bioconductor
#>  vctrs                  0.2.4    2020-03-10 [1] CRAN (R 3.6.0)
#>  withr                  2.4.1    2021-01-26 [1] CRAN (R 3.6.0)
#>  xfun                   0.15     2020-06-21 [1] CRAN (R 3.6.0)
#>  XML                    3.99-0.3 2020-01-20 [1] CRAN (R 3.6.0)
#>  xml2                   1.3.0    2020-04-01 [1] CRAN (R 3.6.0)
#>  XVector                0.26.0   2019-10-29 [1] Bioconductor
#>  yaml                   2.2.1    2020-02-01 [1] CRAN (R 3.6.0)
#>  zlibbioc               1.32.0   2019-10-29 [1] Bioconductor

@Alex Does this solution work when aligning multiple genome tracks using ggbio as in my original example?  Either way, useful info for others. My original question was actually aiming to get feedback on how others make these types of plots rather than seeking a solution to a specific problem in my workflow, so if you have any thoughts on that I'd be interested

@Reilstein this solution requires more work to align multiple genome tracks based on coordinates, so actually does not really answer your initial question.
I was using Gviz before, but I like ggbio more due to the ease of use and how it integrates with ggplot2.


text:
Biopython Phylogenetic Tree replace branch tip labels by sequence logos
Having recently constructed a lot of phylogenetic trees with the module TreeConstruction from Phylo package from Biopython, I've been asked to replace the branch tip labels by the corresponding sequence logos (which I have in the same folder). I thought that it would be more efficient to make a code to generate the logo-trees automatically, as I would have to make a lot of them.
The first idea that I came up with was to see whether the functions used to build the tree had an argument to replace the branch tip labels or to remove them, which I couldn't find. Therefore, I removed the branch tip labels by setting their font size to 0: (following is the code to build the tree)
# Modules to build the tree
from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor
from Bio.Phylo import draw
from Bio import Phylo, AlignIO
import subprocess
import matplotlib
import matplotlib.pyplot as plt

alignment = AlignIO.read('MotifSeqAligned.fasta', 'fasta') # reading the alignment file

calculator = DistanceCalculator('ident')
dm = calculator.get_distance(alignment) # distance matrix

constructor = DistanceTreeConstructor()
tree = constructor.nj(dm) # build with neighbour joining algorithm a tree from dm

Phylo.write(tree, 'TreeToCutOff.nwk', 'newick')

plt.rc('font', size=0)          # controls default text sizes #HERE IS THE SETTING FOR THAT ALLOWS ME TO HIDE THE BRANCH TIP LABELS
plt.rc('axes', titlesize=14)     # fontsize of the axes title
plt.rc('xtick', labelsize=10)    # fontsize of the tick labels
plt.rc('ytick', labelsize=10)    # fontsize of the tick labels
plt.rc('figure', titlesize=18)   # fontsize of the figure title

draw(tree, do_show=False)
plt.savefig("TreeToCutOff.svg", format='svg', dpi=1200)

From this code I could get the tree:

Since I don't know how to get the y coordinates of the branches to add the logos one by one, I built a column of logos with matplotlib, that I intended then to paste on the tree in python. The code to build the column of logos is the following:
#Extract filename from newick
newickFile = open("TreeToCutOff.nwk", 'r').read()
orderedLogos = ["{}.eps".format(i) for i in re.split('(\W)', newickFile) if "Profile" in i]

#Initialize the figure
fig = plt.figure()

#Add each image one after the other in the right order
for i, files in enumerate(orderedLogos):
    img1 = mpimg.imread(files)
    ax1 = fig.add_subplot(len(orderedLogos), 1, 1+i)
    ax1.imshow(img1)
    ax1.set_xticks([])
    ax1.set_yticks([])

# plt.show()
plt.savefig("RowsOfLogos.svg", format='svg', dpi=1200)
plt.clf()
plt.cla()

Having my tree and the column of logos in .svg or .png, I couldn't find any way to stack them properly. My first idea was to use the library svgutils which seemed to be easy to handle, with the following code: (taken from svgutils tutorials)
import svgutils.transform as sg
# Assemble
#create new SVG figure
fig = sg.SVGFigure("14cm", "14cm")

# load figures
fig1 = sg.fromfile('TreeToCutOff.svg')
fig2 = sg.fromfile('RowsOfLogos.svg')

# get the plot objects
plot1 = fig1.getroot()
plot2 = fig2.getroot()
plot2.moveto(280, 100, scale=0.05)

# append plots and labels to figure
fig.append([plot1, plot2])

But the issue with the output was that the background of the column of logos was white and thus, I was pasting a huge white image with a thin column of logos on the tree. And I couldn't find a way to crop the column of logos with svgutils. I tried the module Image from PIL package to build a tree of logos from .png files but couldn't see the tree used as background after pasting the column of logos.
There may be a way to do what I'm aiming for with matplotlib (which would be to stack 2 files .png, and place the logos all the time at the same distance), but I couldn't work it out.
Does anyone know what the best solution is to make a tree of logos (as in the following image which I could only build manually with inkscape) with python libraries, allowing to automate the process without having to adapt the code depending on the number of branches?

Following is a subset of "MotifSeqAligned.fasta" containing the aligned sequences used to build the trees:
>ProfileCluster0.meme
---SSNDTTTCCAGGAAD-
>ProfileCluster1.meme
YBNRD---TTCYYGGAAT-
>ProfileCluster10.meme
---VDKDWTTCTYGGAAT-

With the 3 corresponding logos:

The full "MotifSeqAligned.fasta" and all the logos.eps (as I used them instead of .png which is the format asked by the forum) can be found here.

Thank you for putting so much effort into your question! If I could ask for one more thing though, it would be great if you could give us a minimal, complete and verifiable example, code that we can use directly on our system to replicate your issue so we can play with it and see if we can find a solution.

Thank you for the suggestion, I added the inputs at the end of the topic. Note that the logos are .eps files and here I added them as images (.png). The full inputs are on the link at the end of the topic.

If there's no obvious biopython solution have a look at the ETE Toolkit

To expand on my comment from yesterday. You could do this with the ETE Toolkit (I just copied one logo file rather than converting all 26 to png):
from ete3 import Tree, TreeStyle, faces
def mylayout(node):
    if node.is_leaf():
        logo_face = faces.ImgFace(str.split(node.name, '.')[0] + ".png") # this doesn't seem to work with eps files. You could try other formats
        faces.add_face_to_node(logo_face, node, column=0)
    node.img_style["size"] = 0  # remove blue dots from nodes

t = Tree("tree.nwk", format=3)
ts = TreeStyle()
ts.layout_fn = mylayout
ts.show_leaf_name=False # remove sequence labels
ts.scale =  10000 # rescale branch lengths so they are longer than the width of the logos
t.render("formatted.png", tree_style = ts, h=3000, w=3000) # you may need to fiddle with dimensions and scaling to get the look you want

If you want all of the logos lined up in a column add aligned=True to faces.add_face_to_node

Thanks! That's exactly what I was looking for! I will try with all my logos and let you know if I have further issues.


text:
Error creating indices using STAR
I am trying to index wheat genome using STAR through following command
STAR --runMode genomeGenerate --genomeFastaFiles Triticum_aestivum.TGACv1.dna_sm.toplevel.fa --runThreadN 28

But getting following error,
terminate called after throwing an instance of 'std::bad_alloc'
what():  std::bad_alloc
Aborted (core dumped)

The size of the genome is 13.5 GB and the server capacity is 32 cores with 125 GB RAM, how this error can be resolved?
Memory usage:
I have decreased the number of threads to 8 but again getting the same error as shown below:
> STAR   --runMode genomeGenerate   --runThreadN 8   --genomeDir wheat   --genomeFastaFiles wheat.fa
Jun 04 02:13:31 ..... started STAR run
Jun 04 02:13:31 ... starting to generate Genome files
terminate called after throwing an instance of 'std::bad_alloc'
what():  std::bad_alloc
Aborted (core dumped)

Even I have tried to use STAR without any threads but still the same error:
> STAR   --runMode genomeGenerate   --genomeDir wheat   --genomeFastaFiles wheat.fa
Jun 04 02:19:44 ..... started STAR run
Jun 04 02:19:44 ... starting to generate Genome files
terminate called after throwing an instance of 'std::bad_alloc'
what():  std::bad_alloc
Aborted (core dumped)

May there can be an issue in reference of wheat that I am using, I got this reference from here.
Here are the first ten lines of the reference using head command:
> head wheat.fa
>Pt dna:chromosome chromosome:TGACv1:Pt:1:135835:1 REF
ACAGAAATACCCAATATCTTGTTCTAGCAAGATATTGGGTATTTCTGTCTTTTCTTTCTT
CAAAAATTCTTATATGTTAGCAGAAAAACCTTATCCATTAATAGATGGAACTTCAACAGC
AGCTAAGTCTAGAGGGAAGTTGTGAGCATTACGTTCGTGCATTACTTCCATACCAAGGTT
AGCACGGTTGATGATATCAGCCCAAGTATTAATAACGCGACCTTGACTATCAACTACAGA
TTGGTTGAAATTGAAACCATTTAGGTTGAAAGCCATAGTACTAATACCTAAAGCAGTGAA
CCAGATTCCTACTACAGGCCAAGCAGCCAAGAAGAAGTGTAAAGAACGAGAGTTGTTGAA
ACTAGCATATTGGAAGATTAATCGGCCAAAATAACCATGAGCAGCCACAATATTATAAGT
TTCTTCCTCTTGACCAAATTTGTAACCCTCATTAGCAGATTCATTTTCAGTAGTTTCCCT
GATCAAACTAGAGGTTACCAAGGAACCATGCATAGCACTGAATAGGGAACCGCCGAATAC

How this problem can be resolved?

Keep an eye on the memory usage and see if you're actually exceeding what's available.

The memory does not exceeds 1.6 usage that is normal,

Here is the link to the memory usage

Try removing the --sjdbGTF... stuff from your command as it's likely there's an issue either with the GFF3 file or the fasta file.

Use (much!) fewer threads. Threads speed up processing, but additional threads increase the memory requirement linearly.

@KonradRudolph, why not put your comment as an answer? I would certainly upvote if it was an answer.

@b.nota Fair enough, upon reflection it would probably solve OP’s problem.

In addition to the other answer (to reduce the number of threads), you might also need to add a --genomeChrBinNbits argument. Your reference has more than 5000 scaffolds, so you need to adjust the --genomeChrBinNbits argument to reduce RAM consumption. Found here.
They recommend:

--genomeChrBinNbits} = min(18, log2(GenomeLength/NumberOfReferences))

In your case:
Genome size 17G, number of scaffolds 735945
So --genomeChrBinNbits should be set to log2(17,000,000,000/735,945) = 14.4955, so let's say 14 or 15.

After your command it gave error to increase ram using limitGenomeGenerateRAM, then I increased the ram using limitGenomeGenerateRAM and finally it worked fine, thanks.

@AmmarSabirCheema, I am glad it finally works! Thanks for your feedback.

Try reducing the number of threads.
When multithreading the index creation (and other memory bound tasks), the memory usage increases linearly with the number of threads. If you required a 10 GiB working space work with a single thread, then the requirement for 10 threads would be at 100 GiB.
STAR is particularly memory hungry for the index creation to start with, and your genome is exceptionally large. I’ve had STAR require > 32 GiB of memory for the human genome, with a moderate number of threads. I’m not surprised that it fails with a 4× larger genome and many more threads: it would probably require at least 2× as much RAM as your machine has.

Thanks for your answer but it didn't work.

I have again edited my post to add more details to my problem.


text:
Pathway level analysis of single-cell gene expression
I'm looking for single-cell specific methods to construct (using gene expression data) new features that express pathway "level" or "activity", and then use these for clustering cells.
One example for bulk RNA-seq is PLAGE, implemented in the GSVA R package.
I have found only one method for single-cell RNA-seq that comes close: pagoda in the scde R package. This performs PCAs for each pathway, using its genes.
A very similar question is Classification (supervised learning) of expression data on pathway level, on bulk/microarray data.

Hi, I find a bit unclear your question, could you clarify if you want to find (new?) pathways in scRNA-seq or how to analyze scRNA-seq using pathways? If the later, what is the problem with the methods you mention?

I'm not trying to find new pathways, I want to use pathway- or biological function-level information to classify or understand cell diversity. GSVA is designed for bulk data, so I don't yet know how applicable it is to single cells. scde might be good.

1)Pathways might (most probably do) change for each cell line (and in fact if it is not from a isolation, you don't know if this is a different cell line or not unless you traced it somehow) 2) We don't know enough about pathways and how they work (ie, there is few reliable information available, different databases define different what a pathway is, either by interactions, with kinetic constants an thus by contact (but then in which conditions are these measured?))

Thanks for pointing me to scde. I found that it identifies "known pathways or novel gene sets"

Package Hipathia allows you to compute a value of "activity" for each pathway and each cell, so that the matrix of gene expression is transformed into a matrix of pathway activity. You can use this matrix for further analysis, cell differentiation and clustering. You can also compute functional activity matrices in the same way.

The AUCell R package identifies enriched gene sets within cells. It uses the Area Under the Curve (AUC) to calculate whether a gene set is enriched within the expressed genes for each cell.

Single-Cell_Signature_Explorer paper compared several such methods including: AUCell, Seurat CellCycleScore, and GSVA/ssGSEA.

Check our recent paper where we, not only characterized pathway activities in scRNAseq experiments, but also took advantage of mechanistic models in order to make some in silico prediction of drug effects.


text:
Smallest group size for differential expression in limma (bulk RNA-Seq)
I am reading Smyth et al. (ref. 1). I want to run differential expression analysis on a bulk RNA-Seq dataset in which each group is composed by 2 samples. In the paper previously cited it is written that:

Genes must be expressed in at least one group (or in at least three samples across the entire experiment, where three was chosen as this is the smallest group size) to be kept for downstream analysis.

Is it possible to use limma DE analysis also with groups composed by only 2 samples? NB. It is possible that in this particular dataset the smallest group size is 3.
If not, which alternative should I use?
Update
I have 17 samples. My idea is to test one group vs the rest of groups (all together). So at least the reference group would be composed by > 10 samples. In this case, what can be a feasible analysis for DE?

Smyth, G. K., Law, C. W., Alhamdoosh, M., Su, S. & Ritchie, M. E. RNA-seq analysis is easy as 1-2-3 with limma, Glimma and edgeR. F1000Research 5, 1408 (2016).

The problem with most of the methods is that use the gene's variance for each group, which can't be calculated (reliable) when the sample is <= 2. Also, statistically, it would have extremely low power, so the conclusions couldn't be trusted much.
You could also calculate the "raw" fold change by yourself (ie without the variance estimation and adjusting). If you had one more sample per group you could use DESeq2 which was thought for such cases. But from the paper:

However, if there are two or fewer replicates for a condition, these samples do not contribute to outlier detection, as there are insufficient replicates to determine outlier status.

Which makes me unsure if it will work well for just two samples as it won't help to detect outliers genes.

Thanks. I updated my question. In particular, I was thinking about using your suggestion and just take the raw fold-change. Is it correct to take genes that are considered outliers (e.g. > 3 StD) when comparing the group to the rest of samples?

The problem with such a few number of samples is that it is hard to know what is an outlier and what it is "normal". If you compare 2 vs 15 you still have two problems. 1) you still need to calculate the variance for the group of 2, which has the same problem as initially and 2) the comparison will be meaningful? If you considered them to be different groups is with a reason, if you mix them it will help you to understand what are the differences between them?

You could show the samples and the relevant data about these samples (which groups, do they belong, if they are biological or technical replicates or the differences between the conditions, ...) to help you better.

The dataset I'm using is this one on GEO (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE109125). There's no citation yet and I didn't find in the metadata if they are biological or technical replicates. However, I think they are biological replicates. The group is the first part of sample name (before the '#' sign). The full normalized matrix can be downloaded here: ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE109nnn/GSE109125/suppl/GSE109125%5FGene%5Fcount%5Ftable%2Ecsv%2Egz

This may sound obvious, but can you provide some context in the answer on why the variance cannot be calculated reliable when the sample is <= 2?

@gc5 With less than 2 you don't have variance, and with just 2 it is not good enough, as it is defined by the difference between the points and the mean, and if you add another point you add 33% more data! But that would be a good question for stats.SE, who could give a more reasonable and detailed answer than me :\

While I share Llopis' concern about estimating variance from 2 samples, the statement you quoted is about avoiding false positives from genes that are only expressed in a few samples. It's fairly common to exclude genes that expressed in fewer samples than the smallest group even if the number of samples per group is much higher than 3.

Good point, thanks. I also think the sentence meant that in this specific case.


text:
How can I get a list of papers where a particular mutation has been mentioned?
I have recently tried a few tools (MutationFinder, tmVar, PubTator etc.) that extract mutations from a given text. However, for my work I need something that will take a list of mutations as an input and return a list of pubmed id in which the mutations has been found.
Example input:
mut_list = [Y702F, L483Y]
Example output:
pub_list = [12215417, 18537569]
Can anyone give any idea?
Thanks!

Welcome to the site. Could you elaborate more your question? What did you try to solve this problem?

You could use the pubmed API to query it directly.  Typically I use Biopython for my pubmed queries.  There are some good examples in the cookbook and on other SO posts.
This example from the cookbook link above has the guts of what you need.  You will just have to iterate over your list.  Depending on the number of queries you want to do you might have issues getting throttled by NCBI
In this example, we will query PubMed for all articles having to do with orchids (see section 2.3 for our motivation). We first check how many of such articles there are:

>>> from Bio import Entrez
>>> Entrez.email = "A.N.Other@example.com"     # Always tell NCBI who you are
>>> handle = Entrez.egquery(term="orchid")
>>> record = Entrez.read(handle)
>>> for row in record["eGQueryResult"]: ...
        if row["DbName"]=="pubmed": ...
            print(row["Count"]) 463

Now we use the Bio.Entrez.efetch function to download the PubMed IDs of these 463 articles:

>>> from Bio import Entrez
>>> handle = Entrez.esearch(db="pubmed", term="orchid", retmax=463)
>>> record = Entrez.read(handle)
>>> handle.close()
>>> idlist = record["IdList"]

This returns a Python list containing all of the PubMed IDs of articles related to orchids

A follow-up question: What will be the value of the db parameter if I want to search all databases? I read through Biopython's Entrez.esearch documentation as well as source code but wasn't able to figure it out.


text:
Retrieving gene matrix from internet in GSEA desktop
This error appeared to me, and I have the program connected to the internet. I already installed java 8, but it still doesn't get me the gene set databases.

Did you hit the "Help" button at the bottom?

Yes i did, but it didn't help

Which version of the software are you using? Did you check your connection ?  Why do you need an online "Gene Matrix" dataset? Aren't  you trying to analyze your data?

Have you checked firewall settings on your laptop?
See Firewall / FTP connection issues at:
GSEA wiki - known issues


text:
How to split bedfile into non-overlapping regions and compute aggregation function on duplicate segments
I have a bed-like file in which I have several overlapping regions with associated features. I want to split the overlapping regions into disjoint regions in order to process the file (e.g. with pandas), group by regions and compute an aggregate function.
The input data format is the following (overlapping are marked with *):
chr9    8855925    8986425    -1  *
chr9    8855965    8986325     0  *
chr9    8893857    9022105    -1

I want to split into disjoint intervals, then I want to groupby per-interval and compute an aggregation function (e.g. the mean of the last field):
chr9    8855925    8855965    -1
chr9    8855965    8893857  -0.5
chr9    8893857    8986325 -0.67
chr9    8986325    8986425    -1
chr9    8986425    9022105    -1

I was thinking of using bedtools and pandas to aggregate. However, I was unable to find the right tool to do this job in the bedtools/pybedtools documentation.

To make disjoint intervals, you could use BEDOPS bedops --partition, piping to bedmap --mean to get the mean signal over disjoint regions.
Starting with the input bedgraph file, convert it to five-column BED with GNU awk, putting the signal in the fifth column per UCSC convention:
$ awk -vOFS="\t" '{ print $1, $2, $3, ".", $4 }' /tmp/in.bedgraph | sort-bed - > /tmp/in.bed
$ cat in.bed
chr9    8855925 8986425 .   -1
chr9    8855965 8986325 .   0
chr9    8893857 9022105 .   -1

The disjoint set from these intervals would look like this:
$ bedops --partition /tmp/in.bed
chr9    8855925 8855965
chr9    8855965 8893857
chr9    8893857 8986325
chr9    8986325 8986425
chr9    8986425 9022105

We can pipe these intervals to bedmap --mean with the original five-column BED file as a "map" file, calculating the mean signal over the disjoint intervals:
$ bedops --partition /tmp/in.bed | bedmap --echo --mean --delim '\t' - /tmp/in.bed > /tmp/answer.bedgraph
$ cat /tmp/answer.bedgraph
chr9    8855925 8855965 -1.000000
chr9    8855965 8893857 -0.500000
chr9    8893857 8986325 -0.666667
chr9    8986325 8986425 -1.000000
chr9    8986425 9022105 -1.000000

The --mean option calculates the arithmetic mean. Other aggregation functions are available to process the signal or score (fifth) column in map BED files, like --stdev, --median, --mad, etc.
When calculating statistics, the default precision is six digits. The --prec <n> option can be used to specify desired precision in the score result.
See bedmap --help or the online documentation for a full description.
To avoid making and cleaning up intermediate files, bash process substitutions could be used to generate a one-liner:
$ bedops --partition <(sort-bed in.bedgraph | awk -vOFS="\t" '{ print $1, $2, $3, ".", $4 }') | bedmap --echo --mean --delim '\t' - <(sort-bed in.bedgraph | awk -vOFS="\t" '{ print $1, $2, $3, ".", $4 }') > answer.bedgraph

This can be useful for cleaner pipelines or with processing very large, whole-genome scale inputs on systems with limited disk space.

Thanks for the detailed answer. I got confused with the numbers, I should have used very simple segment start and end. Your answer has the expected solution. I'll update my question.

On macOS, use gawk instead of awk since the -v option is invalid. With homebrew brew install gawk.


text:
Raw vs Filtered in the output of cellranger count
After running cellranger count I got two relevant for further analysis folders: filtered_gene_bc_matrices and raw_gene_bc_matrices. What is the difference between them? What does cellranger filter out, and what should we preferably use for further analysis?

Mostly people use the filtered gene bc matrix. Take a look at the guided tutorial for the r package Seurat if you're new to 10x and /or single cell (https://satijalab.org/seurat/pbmc3k_tutorial.html)

I'm unsure whether this is the answer you are looking for, but when looking into 10X cellranger documentation for the Matrices Output:

Unfiltered gene-barcode matrices:
  Contains every barcode from fixed list of known-good barcode sequences. This includes background and non-cellular barcodes.
Filtered gene-barcode matrices:
  Contains only detected cellular barcodes.

Source


text:
How to apply function on bedfile based on a reference bedfile
I have two bed-like files:
File1
chr1     1     10    .    1    Sample1
chr1    20    180    .    1    Sample1
chr1     1      5    .    1    Sample2
chr1    70    100    .    1    Sample3

File2 (reference)
chr1     1     5     .    0
chr1     5    20     .    1
chr1    50   120     .   -1

I'd like to find a way to apply a function (in my case a subtraction) between the first file and the second, using the second as a reference, and producing a results file containing also the indication of the sample name.
For example, in my case the results file would be:
chr1     1      5    .    1    Sample1
chr1     5     10    .    0    Sample1
chr1    20     50    .    1    Sample1
chr1    50    120    .    2    Sample1
chr1   120    180    .    1    Sample1
chr1     1      5    .    1    Sample2
chr1    70    100    .    2    Sample3

I was thinking about implementing it in pandas, but it would be long to do and probably not so efficient. I looked for bedmap but it seems to me that this operation can be carried on (with the sum operation) only within a single file during aggregation, not between two different bed files. So maybe it is not the right tool to use in this case.
Also, it would be important in the final file to retain the sample name.

You could munge your first file into a BED5 file:
$ awk -vOFS="\t" '{ print $1, $2, $3, $6, $5 }' /tmp/f1.bed | sort-bed - > /tmp/f1.bed5

Sort the second file, if necessary:
$ sort-bed /tmp/f2.unknownSort.bed5 > /tmp/f2.bed5

Take the union of the two sorted files:
$ bedops -u /tmp/f1.bed5 /tmp/f2.bed5 > /tmp/union.bed

Here's what the union looks like:
chr1    1       5       .       0
chr1    1       5       Sample2 1
chr1    1       10      Sample1 1
chr1    5       20      .       1
chr1    20      180     Sample1 1
chr1    50      120     .       -1
chr1    70      100     Sample3 1

We have one name scheme for reference element IDs (.) and a second name scheme for sample IDs (Sample*).
We will exploit this difference in name schemes when we apply operations in a certain order. For this example, we will subtract a reference element score from a sample element score.
We can partition the union, mapping the disjoint elements back to the unioned set:
$ bedops --partition /tmp/union.bed | bedmap --echo --echo-map-id --echo-map-score --delim '\t' - /tmp/union.bed
chr1    1   5   .;Sample2;Sample1   0.000000;1.000000;1.000000
chr1    5   10  Sample1;.   1.000000;1.000000
chr1    10  20  .   1.000000
chr1    20  50  Sample1 1.000000
chr1    50  70  Sample1;.   1.000000;-1.000000
chr1    70  100 Sample1;.;Sample3   1.000000;-1.000000;1.000000
chr1    100 120 Sample1;.   1.000000;-1.000000
chr1    120 180 Sample1 1.000000

Each of the fourth and fifth columns are semi-colon-delimited strings. Each element is ordered by the order in which it is discovered in a bedmap operation.
Both columns have the same ordering or same one-to-one correspondence of ID to signal score. This is very useful, because we can turn these into a pair of ordered arrays and apply some logic, depending on the answer we want back.
For example, here is some awk logic to parse those two strings and subtract a reference element score from a sample element score:
$ bedops --partition /tmp/union.bed \
    | bedmap --echo --echo-map-id --echo-map-score --delim '\t' - /tmp/union.bed \
    | awk -vFS="\t" -vOFS="\t" \
        '{ \
             i=split($4,a,";"); \
             j=split($5,b,";"); \
             if ((i==1) && (a[i]==".")) { pass; } \
             else { \
                 for (k=1; k<=i; k++) { \
                     if (a[k]==".") { refIdx=k; } \
                 } \
                 for (k=1; k<=i; k++) { \
                     if (k!=refIdx) { print $0, ".", b[k]-b[refIdx], a[k]; } \
                 } \
             } \
         }' \
    | cut -f1-3,6-8 \
    | sort-bed - \
    > answer.bed

The file answer.bed will look like this:
chr1    1   5   .   1   Sample1
chr1    1   5   .   1   Sample2
chr1    5   10  .   0   Sample1
chr1    20  50  .   1   Sample1
chr1    50  70  .   2   Sample1
chr1    70  100 .   2   Sample1
chr1    70  100 .   2   Sample3
chr1    100 120 .   2   Sample1
chr1    120 180 .   1   Sample1

This partitioning of intervals is different from the one in your question. Perhaps this is due to how you partitioned your set. In any case, perhaps this would give you an idea of how to apply custom operations over a disjoint partitioning.
The main idea is that this part of the awk script finds the score associated with the reference file ("File2"):
for (k=1; k<=i; k++) {
    if (a[k]==".") { refIdx=k; }
}

Any index other than refIdx corresponds to a Sample* element. That means we can do operations between Sample* scores and the specific reference score.
We also leave out mappings where a reference score does not associate with any Sample* element:
if ((i==1) && (a[i]==".")) { pass; }

Thus, all the operations are always between a pairing of a sample and the mapped reference element.
The cut and sort-bed statements at the end just clean up the result so that you get something that looks like what you are expecting at the end, which is also sorted correctly and ready to use for downstream set operations.


text:
RunUMAP in Seurat not working: module 'umap' has no attribute 'UMAP'
I am trying to run UMAP in the following way:
RunUMAP(seurat_object, dims.use = 1:15, reduction.use = "pca", reduction.name = "umap",
    reduction.key = "UMAP", n_neighbors = 30L, min_dist = 0.3, metric = "correlation")

Or just this:
RunUMAP(seurat_object, dims.use = 1:15, reduction.use = "pca")

It is giving me the error:

Error in py_get_attr_impl(x, name, silent) :
    AttributeError: module 'umap' has no attribute 'UMAP'

I already installed umap-learn successfully in this way: conda install -c conda-forge umap-learn and checked that Rstudio is using the same python version as is provided by default in the terminal: 3.6.5

You should post this on their github page as an issue, looks more like a software use or bug question rather than a bioinformatics question.

This is resolved by installing the more recent version of the code from GitHub.
git clone https://github.com/lmcinnes/umap
cd umap
pip install --user -r requirements.txt
python setup.py install --user

As Peter said, this kind of question is best directed to the Issues page on their GitHub repository.


text:
Why Ti/Tv ratio?
I'm interested in the transition/transversion (Ti/Tv) ratio:

In substitution mutations, transitions are defined as the interchange of the purine-based A↔G or pryimidine-based C↔T. Transversions are defined as the interchange between two-ring purine nucleobases and one-ring pyrimidine bases.

(Wang et al, Bioinformatics, 31 (3), 2015)

How exactly does this ratio imply false positives? Too high ==> high false positive rates? Or too low?
Why is the expected value for random substitutions for the Ti/Tiv ratio 0.5?
If the ratio is expected to be 2.10 for WGS, but I get 3.00 what does that mean? What if I get 1.00?

I know nothing about this ratio (I had to search to know that it is about transitions to the number of transversions), but this link might help others.

I more like to use "ts/tv" for transition-to-transversion ratio. This abbreviation had been used in phylogenetics. When NGS came along, some important developers started to use "ti/tv", but I am still used to the old convention.

Why is the expected value for random substitutions for the Ti/Tv ratio 0.5?

There are six types of base changes. Two of them are transitions: A<->G and C<->T and the other four types are transversions. If everything were random, you would expect to see twice as many transversions – ts:tv=2:4=0.5.

If the ratio is expected to be 2.10 for WGS

The expected ratio is not "2.10 for WGS". It is 2–2.10 for human across the whole genome. You see this number when you align the human genome to chimpanzee or when you focus on an accurate subset of human variant calls. However, in other species, the expected ts/tv may be very different. Also, this number is correlated with GC content. You get a higher ts/tv in high-GC regions, or in coding regions which tend to have higher GC, too. Partly as a result, it is hard to say what is expected ts/tv in exact.

but I get 3.00 what does that mean? What if I get 1.00?

If you get 3.00, your callset is highly biased. If you get 1.00, your callset has a high error rate. Suppose the exact ts/tv is $\beta$ and you observe ts/tv $\beta'\le\beta$, you can work out the fraction of wrong calls to be (assuming random errors have ts/tv=0.5)
$$
\frac{3(\beta-\beta')}{(1+\beta')(2\beta-1)}
$$
This is of course an approximate because $\beta$ is not accurate in the first place and because errors are often not random, so their ts/tv is not really 0.5.

How exactly does this ratio imply false positives? Too high ==> high false positive rates? Or too low?

Too low ==> high false positive rate; too high ==> bias. In practice, you rarely see "too high" ts/tv.

I haven't done variant calling myself, but 5 seconds of googling lead me to this:

This metric is the ratio of transition (Ti) to transversion (Tv) SNPs. If the distribution of transition and transversion mutations were random (i.e. without any biological influence) we would expect a ratio of 0.5. This is simply due to the fact that there are twice as many possible transversion mutations than there are transitions. However, in the biological context, it is very common to see a methylated cytosine undergo deamination to become thymine. As this is a transition mutation, it has been shown to increase the expected random ratio from 0.5 to ~2.01. Furthermore, CpG islands, usually found in primer regions, have higher concentrations of methylcytosines. By including these regions, whole exome sequencing shows an even stronger lean towards transition mutations, with an expected ratio of 3.0-3.3. A significant deviation from the expected values could indicate artifactual variants causing bias. If your TiTv Ratio is too low, your callset likely has more false positives.

Since methylated cytosine deamination is a fairly universal phenomenon, if you're seeing a ratio that deviates substantially from this expectation, you've either discovered some exciting new biology, or there are issues with your variant calling pipeline.

The link is broken. Possibly the unbroken link is: https://gatk.broadinstitute.org/hc/en-us/articles/360035531572-Evaluating-the-quality-of-a-germline-short-variant-callset


text:
Can I run STAR without an annotation file?
I wish to use Rascaf to scaffold a fragmented draft genome.
For this, I need to provide a BAM file of aligned RNA-seq reads and the draft genome.
So, I indexed the draft genome with STAR like this:
STAR --runMode genomeGenerate --genomeDir output/index --genomeFastaFiles draft.fa

Then, I tried to aligned the reads like this:
STAR --genomeDir output/index --readFilesIn reads.fastq.gz --readFilesCommand gunzip -c --outFileNamePrefix output/alignment --quantMode GeneCounts --outSAMunmapped None --outSAMtype BAM SortedByCoordinate --outSAMattrRGline ID:RG1 CN:yy \"DS: z z z\" SM:sample

The latter commands outputs the following error:
Transcriptome.cpp:51:Transcriptome: exiting because of *INPUT FILE* error: could not open input file exonGeTrInfo.tab
Solution: check that the file exists and you have read permission for this file
          SOLUTION: utilize --sjdbGTFfile /path/to/annotantions.gtf option at the genome generation step or mapping step

It suggests I should use a GTF file. But this is a draft genome of an individual for which no GTF is available. I could try to adapt a GTF from a close relative but for the scaffolding, it seems superfluous. Is there a way I can map with STAR without the GTF?
These commands were run on a GNU/Linux machine.

I don't think you can use the --quantMode GeneCounts option with no annotations. I think the error is trying to look for an exon file generated from the annotations to do the quantitation on. Remove that and I think it should work, as the manual specifically states that annotations are optional but highly recommended.


text:
How to extract the protein fasta file from a genbank file?
Let's say I have a genbank file, e.g. the one downloaded from here, that contains entries as
 gene            190..255
                 /gene="thrL"
                 /locus_tag="b0001"
                 /gene_synonym="ECK0001; JW4367"
                 /db_xref="EcoGene:EG11277"
                 /db_xref="GeneID:944742"
 CDS             190..255
                 /gene="thrL"
                 /locus_tag="b0001"
                 /gene_synonym="ECK0001; JW4367"
                 /function="leader; Amino acid biosynthesis: Threonine"
                 /GO_process="GO:0009088 - threonine biosynthetic process"
                 /codon_start=1
                 /transl_table=11
                 /product="thr operon leader peptide"
                 /protein_id="NP_414542.1"
                 /db_xref="ASAP:ABE-0000006"
                 /db_xref="UniProtKB/Swiss-Prot:P0AD86"
                 /db_xref="EcoGene:EG11277"
                 /db_xref="GeneID:944742"
                 /translation="MKRISTTITTTITITTGNGAG"

and I want to extract the potein fasta from it, how would I do this?
I found a code sample here. When I then run:
from Bio import SeqIO

gbk_filename = 'NC_000913.3.gb'
faa_filename = 'NC_000913.3_warwick.fasta'
input_handle  = open(gbk_filename, "r")
output_handle = open(faa_filename, "w")

for seq_record in SeqIO.parse(input_handle, "genbank") :
    print "Dealing with GenBank record %s" % seq_record.id
    for seq_feature in seq_record.features :
        if seq_feature.type=="CDS" :
            assert len(seq_feature.qualifiers['translation'])==1
            output_handle.write(">%s from %s\n%s\n" % (
                   seq_feature.qualifiers['locus_tag'][0],
                   seq_record.name,
                   seq_feature.qualifiers['translation'][0]))

output_handle.close()
input_handle.close()

it results in
     31     for seq_feature in seq_record.features :
     32         if seq_feature.type=="CDS" :
---> 33             assert len(seq_feature.qualifiers['translation'])==1
     34             output_handle.write(">%s from %s\n%s\n" % (
     35                    seq_feature.qualifiers['locus_tag'][0],

KeyError: 'translation'

How could this be fixed?

One can get it to work by using SeqIO.InsdcIO.GenBankCdsFeatureIterator:
from Bio import SeqIO

file_name = 'NC_000913.3.gb'

# stores all the CDS entries
all_entries = []

with open(file_name, 'r') as GBFile:

    GBcds = SeqIO.InsdcIO.GenBankCdsFeatureIterator(GBFile)

    for cds in GBcds:
        if cds.seq is not None:
            cds.id = cds.name
            cds.description = ''
            all_entries.append(cds)

# write file
SeqIO.write(all_entries, '{}.fasta'.format(file_name[:-3]), 'fasta')

This will then give the desired output (just the first two entries):
>b0001
MKRISTTITTTITITTGNGAG
>b0002
MRVLKFGGTSVANAERFLRVADILESNARQGQVATVLSAPAKITNHLVAMIEKTISGQDA
LPNISDAERIFAELLTGLAAAQPGFPLAQLKTFVDQEFAQIKHVLHGISLLGQCPDSINA
ALICRGEKMSIAIMAGVLEARGHNVTVIDPVEKLLAVGHYLESTVDIAESTRRIAASRIP
ADHMVLMAGFTAGNEKGELVVLGRNGSDYSAAVLAACLRADCCEIWTDVDGVYTCDPRQV
PDARLLKSMSYQEAMELSYFGAKVLHPRTITPIAQFQIPCLIKNTGNPQAPGTLIGASRD
EDELPVKGISNLNNMAMFSVSGPGMKGMVGMAARVFAAMSRARISVVLITQSSSEYSISF
CVPQSDCVRAERAMQEEFYLELKEGLLEPLAVTERLAIISVVGDGMRTLRGISAKFFAAL
ARANINIVAIAQGSSERSISVVVNNDDATTGVRVTHQMLFNTDQVIEVFVIGVGGVGGAL
LEQLKRQQSWLKNKHIDLRVCGVANSKALLTNVHGLNLENWQEELAQAKEPFNLGRLIRL
VKEYHLLNPVIVDCTSSQAVADQYADFLREGFHVVTPNKKANTSSMDYYHQLRYAAEKSR
RKFLYDTNVGAGLPVIENLQNLLNAGDELMKFSGILSGSLSYIFGKLDEGMSFSEATTLA
REMGYTEPDPRDDLSGMDVARKLLILARETGRELELADIEIEPVLPAEFNAEGDVAAFMA
NLSQLDDLFAARVAKARDEGKVLRYVGNIDEDGVCRVKIAEVDGNDPLFKVKNGENALAF
YSHYYQPLPLVLRGYGAGNDVTAAGVFADLLRTLSWKLGV

Is there any way to include the proteins names annotated at the gbk file as the sequence description?


text:
Specific cell type identification in Single Cell Sequencing
In order to define which cell is of which type we need to identify a set of rules, for instance neurons should express one of the following: Thy1, Rbfox3, MAP2, Camk2b, Gad1,Cck, Reln, and should not express any of the following: cd45, Tmem119, CD11b, ... and others. Maybe Rbfox3 should always be expressed in any neurons. So, these rules need to be compiled manually as I understand to identify each cell type, from the literature. I would assume that many people already faced this issue, and maybe came up with software to do that? Are there any software where we would be able to supply a matrix of single cell expression data and have the types assigned with certain degree of certainty?
I know that maybe we can also use heatmaps for that, like:
cell type specific genes - heatmap using rank based approach
But again, this seems too manual work. Besides, I do not think that we can assign different neuronal types based on just a heatmap.

I do not know of such software.
However, I believe this effort is a bit misdirected. The purpose of single-cell sequencing is to get a better understanding of cells; their heterogeneity and functional diversity or developmental / biological processes such as differentiation, using a higher "resolution" method. In other words, if we had the methods you are asking for, there would be no need for the single-cell experiment.
Running your data through some pipeline utilizing previous knowledge has the danger of forcing ideas on the data, rather than seeing what the data tells you: it would be better to try and understand and explain the biologically relevant heterogeneity and diversity in your cells, together with a critical comparison of cell type characteristics (expression of genes, pathways) with knowledge described in literature.
This question is also very similar to the one you previously asked and was answered.
Also good to recognize that your approach has a lot of implicit assumptions, such as: cells can be clearly categorized as +/- expression for each gene, or it ignores the systems level (networks and pathways), thus again forcing concepts on the data.

This said, one method would be to construct a table listing cell type markers as described in the linked answer above, then write a script that determines a cutoff value for expression of genes in your data (see this), then ranks cell types for each cell. For example, you can measure the number or proportion of uniquely genes expressed (although that's very simplistic).
Another option, if you are familiar with machine learning, is to train a classifier on an annotated dataset, and then use that on new data.
Also see a convenience function below, which requires a table of marker genes for each cell cluster (i.e. (the Seurat::FindAllMarkers() output), and a reference df listing genes and corresponding cell types in its HGNC_symbol and Cell_type columns, and returns the table with listing clusters, their marker genes and corresponding cell types.
getCelltypes <- function(markers, reference) {
  marker.celltype <- markers
  marker.celltype$Cell_type <- marker.celltype$gene
  marker.celltype$Cell_type <- with(reference, celltype[match(marker.celltype$Cell_type, gene)])
  return(marker.celltype)
}

Many variants of this reference table and function can be created, feel free to modify it.
In the future, I think we can expect such classifier software that uses reference data from the Human Protein Atlas, HCA and similar projects.

SingleR, a novel computational method for unbiased cell type recognition of scRNA-seq:
https://github.com/dviraran/SingleR

Welcome to the site.  Could you expand your answer to include how to use this tool to answer the question?

singleCellNet is a computational method using supervised machine learning for quantitatively cell type annotation. singleCellNet also enables cross-platform and cross-species comparison. It is also actively-supported.
Here is the github https://github.com/pcahan1/singleCellNet
, and the vignette https://pcahan1.github.io/singleCellNet/

We have recently published a web-based tool called CIPR (cluster identity predictor) that helps with cluster annotations. You can select one of 7 reference datasets (some contains neurons) or upload a custom reference gene expression file for predicting cluster identity. You can read the manuscript for more information about how the algorithm works, but in a nutshell, it compares differentially expressed genes or global expression patterns within your data to known reference datasets and creates quick visual outputs. I find it helpful for interactive analyses where I try different clustering parameters and such.
This article is a good resource that compares different cluster calling algorithms.


text:
How to create a nexus file for BEAST?
I am learning to use BEAUTI and BEAST software. I followed their tutorial and was able to run example data without any problem. Now, I want to do phylodynamic analysis on my data. I have downloaded 70 picornavirus sequences from the gene bank. My questions are :

What tool do I have to use to create Nexus file that is used by BEAUTI ?

How do I enter the information about the date (year ) the virus was isolated and location information which is required for rates and dates (evolutionary) analysis in BEAST, spatial analysis by SpreaD3 ?

Notes on SpreaD3 is a update on SPREAD by Bielejec et al. SPREAD: Spatial Phylogenetic Reconstruction of Evolutionary Dynamics. Bioinformatics, 2011, 27:2910-2912. This uses spatial Bayesian reconstruction based on the phylogenetic output of Beast

Did you ever find a method of solving this? I have the same issue and seems like the internet is silent on methods and/or too complicated for a beginner.

Use Mega X. Upload fasta files, do Muscle or ClustalW alignment and save output as Nexus.

Its easy,

Download the data from NCBI in fasta format (its an option in one of the pull down headers).
Load it into the latest version of MEGA, version 7(?), perform the alignment therein.
Export from MEGA as Nexus format and upload it into Beast

You could of course perform the alignment in Clustalx, Clustalo or Muscle or T-Coffee, (hmmm I could write an alignment program just to name it T-4two (sorry British joke)). If the output option of "nexus" was not present in the alignment program of your choice, then again upload the alignment into MEGA and export it as nexus.
Anyway, I could go into specific detail about the various points and clicks therein, but I think its easy to figure out. Most phylogeny programs will do an automatic fasta->nexus conversion, but anyway.
Its easy, Beast is not trivial to use properly however and ensuring convergence is not at all easy, in addition the multiple simulation a priori models, whilst powerful, give different answers. To try and explain this a different way, it is easy to get an answer in Beast, but difficult to get a robust answer from it and thats the bit I would give serious consideration about. It does amaze how widespread the program is (I guess its the nice GUI), given how complicated it is to obtain robust result.

Just to add you could try the online tool here for the fasta to nexus conversion:  . I've not tried it, its just an after thought. If you do Python, Biopython will do it, but its just as easy to code.

Thanks ! I did using Mega and it works.