text:
SARS-CoV-2 sequence used in the AlphaSeq Antibody Datasets to predict binding affinity
Currently we are building a sequence based deep-learning model to predict binding affinity between antibody and antigen. For this we are training a sequence based model with AlphaSeq Antibody dataset (https://github.com/mit-ll/AlphaSeq_Antibody_Dataset). But in this dataset they haven't provided the sequence for SARS-CoV-2 virus.
So I read the paper and there they have mentioned they have tested against peptide in the HR2 region of the SARS-CoV-2 spike protein and they have given the sequence as "PDVDLGDISGINAS". (https://www.nature.com/articles/s41597-022-01779-4#Sec3)
But this sequence is very small. I'd like a larger protein sequence that includes non-epitope portions so that the prediction is more meaningful (e.g. predicts binding to the actual epitope when provided a sequence which includes sequence beyond the epitope). I searched online and only found a covid genome sequence which isn't very helpful.
Searching the RCSB PDB website for Spike protein HR2 gave two SARS-CoV-2 results (RCSB PDB website).
One of them had sequence like the following:
>8CZI_1|Chains A, B, C|Scaffolded Spike protein S2' HR1|Nostoc punctiforme (strain ATCC 29133 / PCC 73102) (63737)
MSHHHHHHGSQTLLRNFGNVYDNPVLLDRSVTAPVTEGFNVVLASFQALYLQYQKHHFVVEGSEFYSLHEFFNESYNQVQDHIHEIGERLDGLGGVPVATFSKLAELTCFEQESEGVYSSRQMVENDLAAEQAIIGVIRRQAAQAESLGDRGTRYLYEKILLKTEERAYHLSHFLAKDSLTLGFAYENQKLIANQFNSAIGKIQDSLSSTASALGKLQDVVNQNAQALNTLVKQLSSNFGAISSVLNDILSRLDKVE

>8CZI_2|Chains D, E, F|Spike protein S2' HR2|Severe acute respiratory syndrome coronavirus 2 (2697049)
KNHTSPDVDLGDISGINASVVNIQKEIDRLNEVAKNLNESLIDLQ

So I am not sure which one is the correct sequence for our task. Can someone provide some insights regarding this matter?
Also why are there so many sequences popping up when I search for SARS-CoV-2? Is there a version that can be considered as the root from which other versions are derived?

Can you please replace the image with the actual sequence? The image is hard to read, impossible to copy from and just makes understanding your question harder. You can use the {} button to format it as code.

Is there a reason you don't just use the peptide sequence that they indicate as an antigen? if you're trying to replicate their results then why do you need a different sequence?  Note that your RefSeq COVID link is to the genome (DNA) not protein sequence.

@MaximilianPress Yeah that's a valid question. Actually I wanted to input a bigger sequence to our model rather than the specific part of the whole sequence to which it actually binds, so that our machine learning model learns to guess the correct region from the sequence as well as the correct binding affinity.

@MaximilianPress thanks for indicating that's a DNA

Thanks for explanation, I've edited your question to include this.

PDB has a whole page talking about various slices of the covid spike protein. I'd strongly recommend reading that for a better understanding of why people choose different pieces to represent. In that page they include several different solved structures of the protein that may be helpful for your benchmarking.
In terms of the "canonical" sequence, you could look at a reference genome sequence. Here for example is the EMBL reference genome for covid in GenBank format, which includes the protein sequences. I believe in that the spike protein has ID QHD43416.
You can additionally look at NCBI, which has a slightly different version of the reference (though they should be functionally identical, with just slightly different identifiers and annotations).
Assuming that you want the ancestral initial isolate of covid, you can find a similar GenBank file here from NCBI.


text:
IGV GenArk assembly error
I am simply trying to follow the documentation here: https://igv.org/doc/igvjs/#Reference-Genome/ and test UCSC GenArk assembly GCA_000002305.1 with the following code (tracks removed to avoid complexity)
import React, { createRef, Component } from "react";

const igvStyle = {
  fontFamily:
    "open sans,helveticaneue,helvetica neue,Helvetica,Arial,sans-serif",
  paddingTop: "60px",
  margin: "5px",
};

class IGVViewer extends Component {
  constructor(props) {
    super(props);
    this.container = createRef(null);
  }

  componentDidMount() {
    // Load igv.js dynamically from a CDN
    const script = document.createElement("script");
    script.src = "https://cdn.jsdelivr.net/npm/igv@2.9.4/dist/igv.min.js";
    script.async = true;
    script.onload = () => {
      const { geneData } = this.props;

      // Extract chromosome, start, and end from geneData
      let chromosome = geneData?.chromosomes?.[0] || "NC_051843.1"; // Default chromosome for testing
      const start =
        geneData?.annotations?.[0]?.genomic_locations?.[0]?.genomic_range
          ?.begin || 100000;
      const end =
        geneData?.annotations?.[0]?.genomic_locations?.[0]?.genomic_range
          ?.end || 200000;

      // IGV Options, using GenArk genome reference
      const igvOptions = {
        genome: "GCA_000002305.1", // Replace with any valid GenArk genome ID
        locus: `${chromosome}:${start}-${end}`, // Set locus dynamically from geneData
        tracks: [], // No tracks for now, can add later
      };

      // Log the full IGV options object to verify URLs and locus
      console.log("IGV Options:", igvOptions);

      // Initialize IGV browser
      window.igv
        .createBrowser(this.container.current, igvOptions)
        .then((browser) => {
          console.log("IGV Viewer successfully loaded.");
          // Log the initialized IGV browser to verify further
          console.log(browser);
        })
        .catch((error) => {
          console.error("Error loading IGV:", error);
        });
    };
    document.body.appendChild(script);
  }

  render() {
    return

;
  }
}

export default IGVViewer;

but I get Unknown genome id: GCA_000002305.1 and I am not sure why.

Can you say anything about where this code came from? I followed the tutorial link and didn't see an obvious source. One thing I'll say is that your chromosome NC_051843.1 is a dog chromosome and your genome assembly GCA_000002305.1 is horse. I would not expect that gene to be in that assembly, but I honestly have no idea whether it even got far enough for that to be the issue because I am not used to querying IGV in this fashion. Maybe first try your example in the browser and ensure it works?

source: Stack Exchange

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

sort: votes


text:
I have a partially described gene from a plant, which tool may I use for predicting regolatory and UTRs sequences from a single DNA fragment?
I have a partially described gene from a plant with CDS and 3'UTR described but no 5'UTR; if I know the locus and upstream and downstream regions, is there a tool I can use for predicting regulatory DNA and UTRs for a given gene?
My annotation data source is the assembly of Phaseolus vulgaris here.

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

I think you need to tell us whether the genome has been sequenced. If not the nearest species with a sequenced genome.

when you say "partially described" what do you mean? what's the name? Do you have any sources of e.g. GFF3 annotation files, a la M__'s question?

Hi, sorry, this is my first post here.  So, the plant is Phaseolus vulgaris and is sequenced yet. When I say partially described I mean that I don't have (or can't find) regulatory sequences and in the gff3 file there is not a 5' (the mRNA starts with the first codon). https://ftp.ensemblgenomes.org/pub/plants/release-55/gff3/phaseolus_vulgaris/                    Thank you for the support, I'm new to bioinformatics.

updating with notes from comments

Assuming that you know the gene name, you would be able to look it up in a transcriptome assembly presumably, or perhaps NCBI has a full length transcript (search txid3885[Organism:exp] with your gene name).
For nearby regulatory DNA, you could consult an ATAC-seq dataset or other similar, though I doubt that causal links of regulatory DNA to your gene will have been directly assayed.

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


text:
Would duplicate climatic variables impact maxent?
I want to use maxent to make some future distribution predictions. I downloaded future climate data from worldclim. In worldclim's page there are three variables (max temp, min temp and precipitation) as well as a fourth titled "bioclimatic  variables" which mentions that it contains the other three and some other new ones. My question is: if i download all four and use them for my analysis, will maxent produce reliable results? Or would the duplicate variables be a problem?

When you say "maxent" and "distribution", what do you mean exactly? For example there is specific ecological software called "MaxEnt". But that's also a generic name for maximum entropy methods, which are fairly diverse.

Looking at the data source (which I believe to be this), it does not seem to line up with what you describe. There are more than 3 variables in each of the 4 datasets except the monthly values which has 3, but is presumably a bunch of time series. I am not sure that you will be able to meaningfully combine the time series data with e.g. bioclimatic point estimates for each location (?). You might think to yourself, "would it be handy to use a dataset that someone already preprocessed to be useful for biological modeling?"

Ok I understand better now so I'll write an answer.
Look at the bioclimatic data site. You can download them from the main page, but the suggestion to generate them from the time series data using an R package (dismo). The bioclimatic variables are thus various linear combinations of the 3 variable time series input.
Not having run this myself, I speculate, but it doesn't make sense to me that the bioclimatic variables are also time series, they appear to be summaries across the time series. So they're not exactly identical in information content, but they're close enough that I'm not sure it's worth the effort of figuring out how to munge the monthly data into a format that works with the summaries.
Someone already went to a lot of trouble to think about the biologically relevant representation of this data (e.g. what ANUCLIM does). I recommend relying on that work, or at the least looking at their arguments and deciding whether you agree enough to use their summaries.
It doesn't particularly matter what software you use, I'd argue (probably MaxEnt the species distribution package, but you haven't clarified). Most methods are going to be doing similar things with this data. They are going to find the variables with the best explanation and weigth them more heavily. Throwing in a couple extra correlated variables might get you a little margin, but it is unlikely to increase explanatory power.
If you want a formal methodology to compare the two, you could use something like AIC or BIC to compare the models formally with the 2 variable sets. But it hardly seems worth the effort of the data munging.
Assumption:
I am assuming that you are not talking about adding ~identical variables also into your model, e.g. files here, but are instead talking about the less identical monthly data here. I'm pretty sure this is right, but I wanted to clarify.


text:
Phylogenetic tree building from proGenomes database for shotgun metagenomics
For some weeks I'm fighting with an issue about phylogenetic tree building to use in a phyloseq object in order to calculate beta-diversity metrics that takes into account phylogenetic distance.
I’m trying to create a phylogenetic tree from the representative genomes of proGenomes v.2.1 database. However, I do not have a representative genome for each taxa id. In contrast, I have a multifasta files for more than 100 contigs for each taxa ID (file freeze12.contigs.representatives.fasta.gz).
For each taxaID (the first numbers in each header of the multifasta file is NCBI taxon ID), I would like to bin the contigs into the biggest scaffold possible, and then perform an alignment (through MUSCLE as examples) followed by RaXML-NG to construct a Newick tree.
However, I cannot find any software for binning contigs with only the contigs as inputs. MetaBAT needs bam files from shotgun read mappings to the contigs in question (not present in the proGenomes database). I found another software (Scaffold_builder), but it needs a representative reference genome. And in microbiology I have more than 100 complete genomes for each NCBI tax id...
For this reason, I ask here an open question if exists some software that can perform of contigs taking into account only a multifasta file of contigs?

I would not recommend the approach that you describe. I am not sure that I completely understand it. For example, metagenomic binning by e.g. MetaBAT is not a scaffolding operation, which returns a single sequence. Rather, they are methods to separate a pool of contigs into groups with similar sequence characteristics indicating they come from the same genome.
It sounds like you already have a group of contigs from each genome, you just want them to all be part of a single sequence in a FASTA file, as would be suggested by your interest in the reference-guided scaffolding approach scaffold_builder.
I would propose that you instead use a tool like GTDBtk to place each of your genomes (each represented by a multifasta of contigs, which is the expected input) into an existing taxonomic tree, and then use their placements in that tree as a model of the taxonomy of your genomes, rather than trying to infer it de novo. I assure you that it is a high quality phylogeny that no one will criticize, and it avoids the possibility of introducing errors into your phylogeny through problematic scaffolding operations.
I would suggest reading GTDBtk documentation for more information.
A narrower answer to the original question
I realized that there is an answer to the question of binning contigs without any additional information. I do not recommend it as I think it would give you positively wrong answers because you are starting with your genomes already binned, and it would not provide you with scaffolds, which I think is what you actually want. But there is a way to perform this task.
You can compute tetranucleotide frequencies from contigs and use them to cluster arbitrary sets of contigs into genomes or other groupings based on similarities. This has a significant rate of error, but it does more or less work. It's an older approach.

Sorry for the delay Maximilian Press. Thanks for your comments. I will try to use the GTDBTk way. I really appreciate your knowledge. I will try to go ahead and overcome this situation. Many thanks.

Thank you another time @Maximilian Press. Yes, GTDB-Tk it's a robust software but although the merging of archaea and bacteria taxa is still a current research topic, I have succeed with another tool that permits me obtaining an unique tree (including archea and bacteria).
If I work with GTDB-Tk I've obtained one tree for archaea and other for bacteria. Then to use it in phyloseq I must work with two independent phyloseq objects which is nice also but my purpose of building a tree from shotgun is only to permits me to calculate diversities based on phylogeny.
As alternative that works it was the following:
to construct a phylogenetic tree from Kaiju I used the 40 marker proteins that is based on the proGenomes database used by kaiju (https://progenomes.embl.de/data/proGenomes2.1_markerGenes.tar.gz) and then I've used the ETE-Toolkit(http://etetoolkit.org/) software to construct the tree in Newick format that I can integrate into a phyloseq object.
I post here if someone can help my strategy. It has been no easy (unless for a dummie as me) to construct a phylogenetic tree from Kaiju's taxonomic classifier.


text:
Assemble Anaeroplasma species genome from metagenomic PacBio data
I have a fasta file containing reads generated by PacBio HiFi whole genome sequencing of a feces sample from mouse.
I would like to use this dataset to generate an assembled circularized genome for an Anaeroplasma organism known to be in this sample.
I am familiar with Linux and R. I would very much appreciate your advice on the needed steps and the required tools or R packages.

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.

Many tools prefer PB reads in their original BAM format. What kind of PB reads? HiFi?

Yes HiFi reads. I have the reads in fastq and fasta fromat.

Try hifiasm-meta (paper here). Whether it works depends on the abundance and the complexity. You may also try HiCanu.

I would suggest the following workflow, understanding you to have metagenomic HiFi sequencing that you know to include your species of interest at high abundance:

metagenomic assembly of HiFi with metaFlyE
Annotate your resulting long contigs (>500Kbp) using GTDB-tk
Look for your Anaeroplasma among long contigs
Assess whether your genome is closed.
If the genome is closed, you are done.
If the genome is not closed, you need to close it. Assess its completeness using e.g. checkM.
Map your HiFi reads to the genome, see if you can close the remaining gap with e.g. SSpace-LongRead.

I'd suggest using this paper as an example for further workflow details.

Thanks for your answer and the references

@Amroon, you are very welcome. By the way, it's traditional on SE sites to upvote answers that you find helpful (there are now a couple of reasonable answers to your question), and accept the answer that you see as helping you solve the problem. For more info on responding to answers on SE, see here: https://stackoverflow.com/help/someone-answers


text:
the following errors are found during the installation of an open source software, UNCeqR on z shell
I am trying to do somatic mutation detection via the open source software, UNCeqR
I got following error was detected on mac z shell during installation.
BEGIN failed--compilation aborted at test.pl line 11.
not ok 1
make: *** [test_dynamic] Error 2
  CALLAHAN/Math-CDF-0.1.tar.gz
  /usr/bin/make test -- NOT OK
//hint// to see the cpan-testers results for installing this module, try:
  reports CALLAHAN/Math-CDF-0.1.tar.gz

I would like to ask for any advice to proceed.

Did you look at the cpan-testers results as suggested? Is this the entirety of the output? What command did you run?

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

It is a Perl compilation error:
(Perl libraries) Math::CDF-0.1 -- Generate probabilities and quantiles from several
statistical probability functions

conda install -c bioconda perl-math-cdf

cpan INSTALL Math::CDF

This module provides cumulative probabilities [P(X <= x)] and
quantiles [given p, returns x such that P(X <= x) = p] from several
statistical probability distributions.
This package provides a perl interface to DCDFLIB
https://metacpan.org/release/CALLAHAN/Math-CDF-0.1/source/README
make and make test of Perl libraries Math::CDF-0.1 in the installation are not okay

https://drive.google.com/file/d/1Bd8TklsUPMqDRDn7zImGEXI9Nk84vkfU/view?usp=sharing


text:
Problem with sequence alignment
I'm trying to find if the peptides in column A are present within the longer peptides of column B.
Pairwise alignment tools take the data from both columns as one string (i.e., TIIDYTNNHYTNNHLEII... and PFYMNITGYTDAYNTDIIDYIIDYTNNHLEIITLFN...) which leads to a number of peptides from column A left unaligned even though the peptide is present within one of the longer peptides in column B.
I would like to find a method that works for detecting such peptides.

Easy, few details are you stuck with a CSV format or prepared to copy and paste column A to one file and column B to another. Its solvable either way.

I don't think that you want an alignment so much as you want to test for the presence of an exact substring (column A) within larger strings (column B).
There are simple excel utilities that do this, for example the SEARCH() function. There is no need to use alignments.


text:
Counting the co-occurrence of variants A and B in an aligned sequencing read
I need some help getting started on this project.
To simplify we want to be able to quantify the occurrence of 3 variants on each sequencing read in an alignment file as a proxy measurement for recombination.
The possible combinations we have are A, B or AB.
For example, variant A is at position 100 and variant B is at position 300 of the reference and AB would have both.
We have alignment files against a reference sequence, and then I get stuck. I could technically look by eye in tablet or IGV for example but is there a way to extract each read then query by expected base at the reference position then bin these reads into the categories; A, B, AB and unclassified then calculate the total occurrences of each?
Thanks for your time and suggestions!

You'll need to write something to do this, using pysam or a similar package.

Related question: Access base aligned to particular reference position
After an initial effort posted here, I ended up rewriting this and testing it a little, and posted the script here. It still has some drawbacks (handling indels intelligently) but now it more or less works.


text:
Is there any standardized and ready-to-use dataset available online for download regarding protein geometry?
I want to do ML model training and data analysis for secondary structure prediction.
Is there any standardized and ready-to-use dataset available online regarding protein geometry?
For example, the dataset should have common features like: distances, angles, dihedrals, etc.

I'm confused as to what you are looking for - any data with 2 points will have a distance, and three points will have an angle. Do you mean protein data with the distances, angles, and dihedrals labelled? Why not use a library/tool to calculate these measurements on a dataset?

Your question is not verbose enough to be helpful. What is the specific project and/or problem that you are trying to solve with this? General questions and answers in bioinformatics are rarely helpful to anyone; see https://bioinformatics.stackexchange.com/help/on-topic

@gringer, What is the specific project and/or problem that you are trying to solve with this? --- Secondary structure prediction.

@user366312 please edit and give more details. Secondary structure prediction is not that hard, did you mean tertiary structure? And what exact type of data would you want and what kind of predictions would you be making? Have you read through the SwissModel, or AlphaFold papers? Their references? Are you just looking for Pfam? Or maybe RCSB PDB? We can't help if you don't explain.

I'm closing this question because it is a general question that is not sufficiently described to give good / useful answers. Please edit to add more details, or ask this question on a discussion site like Reddit

As suggested in the comments, PDB is the standard database for protein structural information. I would suggest starting there. They even have some degree of geometric information precomputed for you, there is an API / GraphQL interface that can be used to query certain attributes that do seem to include some geometric measurements and structural annotations.
You will likely have to compute some of the measures you are interested in, it seems that your list is somewhat amorphous regarding which geometric measurements you'd like to measure so I can't recommend much specifically. For very basic things like Ramachandran angles there are various resources you could use, though Biopython is probably the easiest.
There used to be a "Protein Geometry Database" that is now defunct, probably because it is not hard to compute this stuff directly on any given PDB file.


text:
Concatenate multiple fasta files into one file for MLST
I have seven fasta files, one per gene, with more than 400 fasta entries per file, like this:
Input:
Gene1.fasta
>1721_1 gene name
ATG
>10_14 gene name
GTT

Gene2.fasta
>1721_1 gene name
TGA
>10_14 gene name
GAT

Output:
>1721_1
ATGTGA
>10_14
GTTGAT

What I want to do is create one fasta (continuously) file containing seven genes per sample. As you can see, the sample name is not the same before the _ in the fasta headers.
Is there any way to do this with awk or other command-line tools? Thank you very much.

cat Gene*.fasta > all_genes.fasta ?

@TimurShtatland Thank you. But it should contain all genes sequences of each sample. Each sample should be separate files.

I would suggest using samtools faidx to pull the sample's gene out of each file, and then removing the header.  It will ignore anything after whitespace so gene name after space should not be problem. You could iterate across sample names in a bash loop. Of course, it requires that every sample have every gene.

okay. I will try

The concat command from seqkit could concatenate sequences with same ID from multiple files.
seqkit concat [file1] [file2] ...

Assuming your files are small enough to allow you to fit one of them into memory, you can do this:
$ awk '{
        if(/^>/){ name=$0 }
        else{
            if(NR==FNR){seq[name]=$0 }
            else{
                print name;
                printf "%s%s\n", seq[name],$0
            }
        }
      }' file1.fa file2.fa
>1721_1 gene name
TGA
>10_14 gene name
GTTGAT

If you cannot fit the files in memory, you can try this:
$ join -o1.1,1.2,2.2 -t $'\t' <(FastaToTbl file1.fa ) <(FastaToTbl file2.fa) | sed 's/\t//2' | TblToFasta
>1721_1 gene name
ATGTGA
>10_14 gene name
GTTGAT

Note that this approach requires the FastaToTbl and TblToFasta scripts I have posted here: https://bioinformatics.stackexchange.com/a/2651/298.
Or, if brevity is more your thing, the same approach but more compact:

I think it can be done using python dictionary. You can split the file using the header and then make the fasta id as your key and the sequence as your value. Then using simple loop and conditions you can concatenate them.

Okay. I will try

Hi and welcome to the site! Please try and post complete answers. This isn't really answering, it is only giving a general direction in which an answer could be found.


text:
Why does randomized motif search work?
I'm looking for intuition for why a randomized motif search works. My current thinking is as follows:
We are selecting many random kmers from our DNA sequences. The chosen kmers will bias the profile matrix to selecting kmers like them.
Given any particular k-mer chosen, there are two possibilities:

We've selected a meaningless kmer, so the profile will bias selecting random (ie high entropy) kmers and the score will stop improving soon.
We've captured a portion of the true motif. Since the motif isn't random (ie doesn't have high entropy) subsequent iterations using this profile matrix will pick motifs closer to the real motif and the score will keep improving.

Since we give this algorithm many tries and at some point it will eventually arrive at possibility 2 and come to a good conclusion.
Is this a correct understanding of why randomized motif search works? Is there a more rigorous way to prove this?

Update
With regard to comment, the proof of effectiveness of Monte Carlo algorithms like this is a little thorny- if you look at the wikipedia page, you will see that a proof of statistical consistency for any Monte Carlo algorithms took a long time:

From 1950 to 1996, all the publications on Sequential Monte Carlo methodologies, including the pruning and resample Monte Carlo methods introduced in computational physics and molecular chemistry, present natural and heuristic-like algorithms applied to different situations without a single proof of their consistency, nor a discussion on the bias of the estimates and on genealogical and ancestral tree based algorithms. The mathematical foundations and the first rigorous analysis of these particle algorithms were written by Pierre Del Moral in 1996.

Here is a possibly more general reference by Del Moral in 1998. The paper is somewhat dense and I am not an expert in the field so I don't think I'll do it justice, but I think that we could probably figure out a similar mathematical proof for random motif search given time and the right problem formulation.
I sense that your concern might be with Monte Carlo algorithms generally rather than random motif search specifically, so hopefully that gets close to your question.
It does look like people have done heuristic analysis of similar algorithms, but they are not formal proofs.
Original answer
I believe that you are more or less correct.
For example, this presentation's walkthrough of the algorithm (slides 35-36) specifically refers to Greedy randomized profile motif searches. By making the greediness of the algorithm explicit, we can see that there is no claim of finding an optimal solution in any particular iteration, thus the need for multiple retries to get somewhere close to approximate solution.
One additional point made in that presentation is that this algorithm works best when coupled to approaches such as Gibbs sampling for targeting the sampled k-mer chain in regions with higher scores (less entropic k-mers).
A different learning resource with similar content (and a little more intuition) is available here.

That's actually the course I'm taking haha.
I guess what I'm wondering is if we can prove this stuff rigorously. It makes intuitive sense and I'm glad I got the intuition right but I do wonder how we know our intuition is correct.

@Moo see update- I think that addresses the Monte Carlo question, but I can't tell if anyone has done a similar proof for this algo specifically.

Thank you for the update. I think you're right that the root of my question is why Monte Carlo algorithms work.


text:
Kallisto psudoaligner quant not giving me expected output
I am trying to run this kallisto (latest) quant and all I get are abundance and json files.
kallisto quant --plaintext -t 4  -o rf-stranded.nob.1K1_L7.0 -i mm9_ref.idx  --rf-stranded 1K1_L7.0 2K1_L7.0

the results files in stderr is:
[quant] fragment length distribution will be estimated from the data
[index] k-mer length: 15
[index] number of targets: 23
[index] number of k-mers: 363,704,628
[index] number of equivalence classes: 1,945,936
[quant] running in paired-end mode
[quant] will process pair 1: 1K1_L1.0
                             2K1_L1.0
[quant] finding pseudoalignments for the reads ... done
[quant] processed 22,897,424 reads, 4,476,549 reads pseudoaligned
[quant] estimated average fragment length: 24.2427
[   em] quantifying the abundances ... done
[   em] the Expectation-Maximization algorithm ran for 52 rounds

I get an abundance file with:
target_id       length  eff_length      est_counts      tpm
chr1    197195432       1.89745e+08     253015  21.3249
chr10   129993255       1.31621e+08     210320  25.5544
chr11   121843856       1.38776e+08     263606  30.3773
chr12   121257530       1.23195e+08     143568  18.6371
chr13   120284312       1.22484e+08     132969  17.3613
chr14   125194864       1.20926e+08     112711  14.9058
chr15   103494974       1.07866e+08     109335  16.2102
chr16   98319150        9.70075e+07     109045  17.9767
chr17   95272651        1.06375e+08     169166  25.4322
chr18   90772031        9.23542e+07     121990  21.1241
chr19   61342430        7.23377e+07     157142  34.7407
chr2    181748087       1.81791e+08     244870  21.5414
chr3    159599783       1.45978e+08     209257  22.9247
chr4    155630120       1.61978e+08     226893  22.4014
chr5    152537259       1.65672e+08     194951  18.8186
chr6    149517037       1.43243e+08     246652  27.5373
chr7    152524553       1.82795e+08     197527  17.2812
chr8    131738871       1.50556e+08     203954  21.6643
chr9    124076172       1.30446e+08     222639  27.295
chrX    166650296       1.40278e+08     176259  20.0942
chrY    15902555        5.76797e+07     123.462 0.0342311
mtDNA   16299   11804.8 737285  998819
rDNA    249368  720854  33272.5 738.159

run_info.json:
        "n_targets": 23,
        "n_bootstraps": 100,
        "n_processed": 22897424,
        "n_pseudoaligned": 4476549,
        "n_unique": 3641646,
        "p_pseudoaligned": 19.6,
        "p_unique": 15.9,
        "kallisto_version": "0.46.1",
        "index_version": 10,
        "start_time": "",
        "call": "kallisto quant --plaintext -t 4 -b 100 -o rf-stranded.1K1_L1.0 -i mm9_ref.idx --bias --rf-stranded --gtf UCSC_genes.gtf --pseudobam --genomebam 1K1_L1.0.gz 2K1_L1.0.gz"
}

When I ran it with boot strapping (-b), -gtf and --pseudobam --genomebam it gave me a bam  file.
44879990 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
8078785 + 0 mapped (18.00% : N/A)
44879990 + 0 paired in sequencing
22439995 + 0 read1
22439995 + 0 read2
8077160 + 0 properly paired (18.00% : N/A)
8077160 + 0 with itself and mate mapped
1625 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

I'm not sure what is going wrong, can someone help?
NOTE I use a 15nt kmer, and this about 1/2 of what is suggested (31 I think), is that a problem?

Thanks I think I know what is going on.  I probably should give a file of transcripts to the index step and not a chromosome file.  This is why it is giving my only chromosome abundances.

That makes sense. I would suggest posting that as an answer to your own question and then accepting it, to help others.

OK but I first will try to verify.

As you surmised, kallisto requires a transcriptome fasta file, not a genome fasta file.


text:
How to know the number of SNP in my samples after running the population module in stacks
I am new to stacks for analyzing RAD-Seq data, I have run all the modules. I will like to know to get the total SNP from my sample after running the population module. Is there a command line to use and inside which file?
command:
populations -P -M (denovo output) (popmap) -t 5 --vcf --fasta-samples --fasta-loci -O (pop-output) -r 0.8 -p 5

Thanks

Update
Based on the newly provided command (which has --vcf), there should be a VCF file with suffix .vcf in the output directory.
Original answer
It's pretty hard to answer without more information. For example, it's helpful to link to the software and the module that you are working from (I had never heard of Stacks). What do you mean by "total SNP"? Also, what is the actual command that you ran? There are a lot of possible options!
But (based only on reading the docs for the command), it looks like there are a variety of formats in which the SNP data can be output:
File output options:
--ordered-export — if data is reference aligned, exports will be ordered; only a single representative of each overlapping site.
--fasta-loci — output locus consensus sequences in FASTA format..
--fasta-samples — output the sequences of the two haplotypes of each (diploid) sample, for each locus, in FASTA format.
--vcf — output SNPs and haplotypes in Variant Call Format (VCF).
--genepop — output results in GenePop format.
--structure — output results in Structure format.
--radpainter — output results in fineRADstructure/RADpainter format.
--plink — output genotypes in PLINK format.
--hzar — output genotypes in Hybrid Zone Analysis using R (HZAR) format.
--phylip — output nucleotides that are fixed-within, and variant among populations in Phylip format for phylogenetic tree construction.
--phylip-var — include variable sites in the phylip output encoded using IUPAC notation.
--phylip-var-all — include all sequence as well as variable sites in the phylip output encoded using IUPAC notation.
--treemix — output SNPs in a format useable for the TreeMix program (Pickrell and Pritchard).

I would suggest examining the command that you used, specifically what you passed as the output format. Here are the example commands from the docs:
populations -P dir -b batch_id [-O dir] [-M popmap] (filters) [--fstats] [-k [--window_size=150000] [--bootstrap [-N 100]]] (output formats)
populations -V vcf -O dir [-M popmap] (filters) [--fstats] [-k [--window_size=150000] [--bootstrap [-N 100]]] (output formats)

Can you track down what you specified as output? That will contain the SNP data / genotypes.
Possibly you will need to rerun the command while specifying a different, preferred format.


text:
Is it possible to go from a digital sequence to an actual nucleotide?
I'm trying to learn about what the state of technology is at presently.   It seems like we clearly can go from nucleotide to digitally stored sequence, but can we transcribe something from the digital sequence?

Are you asking about synthetic biology? Can you give a bit more context? The answer is it depends, e.g. one can trivially synthesise short primers but not a whole human genome

@Chris_Rands, you mean an amount of material required for an insertion is doable, but you can't synthesize a complete gene?  I'm trying to understand the technological state of bio-engineering, and trying to understand the workflow.

Oligonucleotide synthesis is a service offered by many providers such as IDT.

What is interesting is how phosphoramidite synthesis differs from how polymerases do it. It is a solid state synthesis approach that expands from the 3' to the 5' in cycles of reaction and 2'-deprotection to allow a single nucleotide to be added at a time.

To add to @Maximilian Press' answer.
There is quite a lot of research on how to store arbitrary digital data in DNA. Nick Goldman describes his research on the topic in his TEDx talk. The last progress in the field I have noticed was DNA fountain, a tranformation of bites in DNA with maximal efficiency of data storage in DNA (information / nt) and reliable recovery when sequenced afterwards.


text:
Why do Illumina 850k/EPIC arrays ignore CpGs which are "GC" in the forward strand?
CpGs are symmetrical, in that a CG sequence on the forward strand is hybridized to a GC --- and both dinucleotides on each opposing strand are CpGs dinucleotides which can be methylated. Conversely, CpGs can be GC on the forward strand but CG on the reverse strand.
FORWARD -> 5'--CG--3'  [OR]  5'--GC--3' <- FORWARD
REVERSE -> 3'--GC--5'        3'--CG--5' <- REVERSE

The assignment of "forward" and "reverse" strandedness is more or less arbitrary.
Given the above, why does it seem like the Illumina 850k (aka EPIC) array only profiles methylation from CpGs which are CG in the forward strand, while ignoring CpGs which are GC in the forward strand? I would also love to hear if my premises are wrong.
Here is some R code illustrating the issue at hand:
suppressPackageStartupMessages({
  library(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
  library(tidyverse)})
data(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)

IlluminaHumanMethylationEPICanno.ilm10b4.hg19 %>%
  getAnnotation() %>%
  as_tibble() %>%
  count(forward_seq=str_extract(Forward_Sequence, "\\[[ATCG]{2}\\]"))

# Results:
# # A tibble: 3 × 2
#   forward_seq      n
#   <chr>        <int>
# 1 [CA]          2922
# 2 [CG]        862927
# 3 [CT]            10
```

I cross-posted this to Biostars: https://www.biostars.org/p/9545804/

My question was kindly answered by Biostars user Papyrus, who pointed out that a GC in the 5' to 3' direction is a GpC dinucleotide, not a CpG dinucleotide. GpC dinucleotides are then biologically and chemically distinct from CpG dinucleotides, due to the directionality of the DNA polymer.


text:
How to fix 5'UTR annotation with RNA-seq data?
I am studying the transcriptome of Arabidopsis. Interestingly, the 5' UTR of the latest annotation is usually too long.
Here is an example (this gene is right to left).

You can see both EST (orange bars) and the RNA-seq reads cover a much smaller region of the annotated 5' UTR. Unfortunately, many 5'UTR of the genes have the same problem. Is there a way to fix that? I was wondering if any assembly tools can help me to get the correct 5'UTR ranges. Thanks!!!

What evidence do you have that the gene stops earlier? This observation could just be due to selective amplification of the 3' end in the RNASeq reads.

I think it starts late (not stop early) as I mentioned 5'UTR, not 3'UTR.

Do you really think standard RNAseq is perfect enough to judge the start of a transcript with single-base resolution? There's a reason that people have developed methods explicitly for that.

It is not only my data. It is also from many people's data. There is no read at 5'UTR in the old EST data, and in >100 set of RNA seq data. There is just nothing there. I do not understand why those upstream sequences are annotated as 5'UTR, and it could mislead many sequence analysis.

@Devon Ryan, Could you please enlighten me the reasoning (why did they give a long 5'UTR with no evidence?). I am asking this because there are lots of genes behaving like this. I am afraid it will cause a problem in motif finding. Thanks!
I updated the image to include the EST data.

To be frank, you need to ask the people who made the annotation.

So is there any tool that can help me to reannotate the 5'UTR based on RNAseq data? I am asking this because, as I checked, most cases RNAseq data are quite consistent with the EST data. Thanks!

@l0110 Could you [edit] your comments in the question. I think I would get a better idea of what do you want.

@Llopis, I cannot edit my comments after 5 minutes. Is there any part unclear? Please let me know. Thanks!

In the case of light-grown seedling it looks like there actually are reads (very few) that are actually the entire length of the 5' UTR. I think what you're after is "what is the absolute 5' transcription start site of my genes?" Right? If that is the case then you're better off doing a few SMRT IsoSEQ or ONT cDNA sequencing flowcells to get the starts and stops. Keep in mind that when you select for polyA you are going to get RNAs that have undergone more processing and may have lost some 5' UTR.

Genes, particularly transcribed genes in Eukaryotes, can have alternate start and stop sites (and splicing).  Alternate transcriptional start sites may be found in different tissues or at different developmental periods due to alternate transcription factor binding sites being used.  The only reliable way to determine the transcriptional start and stop sites is by read coverage.
Genome annotations are generally the product of bioinformatic pipelines and are mostly predictions based upon computational models.  Different gene-finders may call different boundaries; for example, in eukaryotic genomes, an exon may be assigned to the wrong up/down-stream gene.  Arabadopsis, as model organism, has had it's gene boundaries considerably refined/corrected/validated, but even so, there could still exist errors in the annotation or, more likely, the boundaries indicated are the most common start/stop sites (see point 1 above).

In short, your UTRs are whatever your experimental data (reads) show it to be.  Don't worry too much about it not matching up to the standard annotation boundaries.

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


text:
Combine Fastq by writeFastq is not working properly
Currently I have 4 fastq file for each patient as WES data. For WES pipeline, I want to combine those fastq files. I am using writeFastq from ShortRead package in Bioconductor.But I think this function is not working properly.
The reason of my doubt is the the size the next step is quality check using FastQC tool. This tool is working on individual file but not on combined file and shows the error that ID line doesn't start with "@" . So I want to know what are the other effective tool for combining multiple fastq files without any loss of information.
Let me know if I am doing something wrong. Thanks
[EDIT-1] : The codes for combining multiple fastq files into one is as follows:
reference_fastq = file.path(paste(getwd(),"refs",sep = "/", collapse = " "), fsep = ".fasta")
fastqDir <- file.path(dataDir, fsep = ".fastq")
outfile_temp <- paste(dirname(dataDir),basename(fastqDir),    paste(basename(fastqDir),"combined","fastq", sep = "."),sep = "/", collapse = " ")
fls <- list.files(fastqDir, "fastq$", full=TRUE)
for (fl in fls) {
  fq <- readFastq(fl)
  writeFastq(fq,outfile_temp, mode="a")
}

First four lines of my fastq file (one of 4 indivual file) are as follows :
@D00239:185:C9KELANXX:2:2201:1088:2107 1:N:0:TCCTGAGC
GTACAAGAAGCATAACATCAGCATGTGTTTTGGTGAGAACCTTAGGAAGTTTGCAATCATGACAGAAGGTGAAGG
+
/<BBBFFFBFFFFFFFFFFFBF<FFBF<B<FFFFFFFFFFB/B/<////<<<FF<B<FF<FFFFFFFFF/FFFFF

First four lines of combines fastq file are as follows:
@D00239:185:C9KELANXX:2:2201:1088:2107 1:N:0:TCCTGAGC
GTACAAGAAGCATAACATCAGCATGTGTTTTGGTGAGAACCTTAGGAAGTTTGCAATCATGACAGAAGGTGAAGG
+
/<BBBFFFBFFFFFFFFFFFBF<FFBF<B<FFFFFFFFFFB/B/<////<<<FF<B<FF<FFFFFFFFF/FFFFF

It is pretty clear that combined fastq file is also starting wiht proper identifier but why FastQC tool is saying that identifier '@' is missing. Here is the screen shot of the error message:

Thanks.

Can you please show the first and last lines of the input FASTQ files, and the first and last line of whatever writeFastq has created?

Could you paste the command you used from ShortRead to combine the files? (And the version of the package, just in case there is a bug in that version). If you have the fastq files compressed you can simply cat them to a new .tar.gz file and you'll have them combined

@Lot, please instead of images paste the text in code formatting. Also could you, answer to gringer?

@gringer: I have shown the first 4 lines of input fastq file and output fastq file. What else you want to check..let me know.

@Llopis: Actually I tried to copy the text instead the image but I am working on server and on right click after text selection, I did not find any option like copy. That's why I took screen shot. Sorry for this.

@Llopis: ShortRead package version is 1.28.0, R version 3.4.3 and bioconductor version is 3.2.

You have a 2 year's old version of Bioconductor, you might be interested in upgrading it (just in case some other errors are already solved in the updated version).

You can drastically simplify your effort because FASTQ files are simple text files; therefore, standard UNIX test file tools work:
cat refs/*.fastq > combined.fastq

And this even works with .fastq.gz files — no unzipping required.

Apart from this, your R code contains several errors. For instance, your variable names change throughout, and you’re using the fsep parameter wrongly. In fact it’s unclear to me what the code is supposed to do. At an approximation (and untested), the following should work:
fastq_dir = file.path(getwd(), "refs")
fastq_files = dir(fastq_dir, pattern = "\\.fastq$", full.names = TRUE)
outfile = file.path(dirname(fastq_dir), "combined.fastq")

for (file in fastq_files) {
    fq = readFastq(file)
    writeFastq(fq, outfile, mode = "a")
}

This is equivalent to running the shell command shown above, just in a more convoluted way (parsing and re-generating the FASTQ files).

@Konard : Thanks for pointing out my mistake. As I am beginner so my codes are not as much good. I have changes this. Can you tell me what is the proper use of fsep. I though using file path, fsep is the way to choose only specific type of file (may be I am wrong). And how can I use cat within for loop because it is giving unexpected symbol error by trying cat fq > combine.fastq Thanks again. I'll try my best to improve. keep pointing my mistakes. :)

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


text:
Multiway intersection of 8 lists of genes
I have 8 lists of differentially expressed genes for 8 time points (2, 4, 6, 8, 10, 12, 14, 16 hours), the number of genes in each of these lists are not equal.
For example there are 630 common genes between h2 and h4 time points but I have 1143 genes in h2 list and 768 genes in h4 list , so common genes in h2 vs h4 would be 630/1143 and common genes h4 vs h2 would be 630/768.
That would would be great if I have a code to produce a matrix or plot of all possible pairwise combinations of my 8 time points.

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

EDIT: Sorry, I did not see the R tag and wrote the code in bash. I still left the bash example below and added an R version.
Depends on what you want to do, but there are two common ways in which you could store this. Either as wide format:
      h2  h4  h6 ...
Gene1 1   1   0
Gene2 0   1   1
Gene3 1   1   1
...

Where 1 denotes presence and 0 absence, or as long format:
      timepoint
Gene1 h2
Gene1 h4
Gene2 h4
Gene2 h6
Gene3 h2
Gene3 h4
Gene3 h6
...

The best format probably depends on what you want to do. Since you did not specify a method/language, I took the liberty to write a few lines of bash to go from the gene lists to the matrices.
Bash examples assume your files are named h2, h4 etc and are tab separated with gene names in the first column. R examples assume tables is a list of your data frames and timepoints is a vector c("h2", "h4", ...).
Generate table in wide format:

Concatenate gene names of all tables to a separate file
Extract unique values in this file
Perform outer join between the separate file and each table

R:
all_names <- c()
for (tp in timepoints){
  all_names <- append(all_names, tables[[tp]]$genes)
}
all_names <- unique(all_names)

wide_table <- matrix(nrow=length(all_names), ncol=(1+length(timepoints)))
wide_table <- as.data.frame(wide_table)
colnames(wide_table) <- append("genes", timepoints)

for (tp in timepoints){
 tmp <- tables[[tp]]
 tmp$tp <- 1
 tmp <- merge(wide_table, tmp, by="genes", all=T)
 wide_table$tp <- tmp$tp
}

Bash:
Note: you will might need to add headers to this table
for timefile in $yourfiles;do
  # make a file with all gene names
  awk -v tp=$timefile '{print $1}' >> names_file
done

sort names_file | uniq > tmp && mv tmp names_file

for timefile in $yourfiles;do
  join -a1 -a2  -j 1 -o 0,1.2,2.2 -e 0 $names_file <(awk '{print$1,1}' $timefile | sort)
done

Generate table in long format:

For each table, add a column with timepoint
Concatenate all tables

R:
df <- data.frame()
for (tp in timepoints){
  tables[[tp]][[timepoint]] <- tp
  df <- rbind(df, tables[[tp]] )
}

Bash:
   for timefile in $yourfiles;do
      awk -v tp=$timefile '{print $1,timefile}' >> long_table
    done

EDIT2: To follow up on this, here is a simple way (using the package dplyr) to count the number of common genes between all possible combination of time points. Assuming you have your data in long format:
df <- data.frame(genes=c("G1","G1","G2","G3","G3"), tp=c("h2","h4","h2","h4","h6"))

  genes tp
1    G1 h2
2    G1 h4
3    G2 h2
4    G3 h4
5    G3 h6

You can get the co-occurrences using :
library(dplyr)
# Make timepoint ordinal factor
df <- mutate(df, tp = factor(tp, ordered=T))
df <- df %>%
  # Compute all pairwise co-occurrences
  full_join(df, by="genes") %>%
  group_by(tp.x, tp.y) %>%
  summarise(N_genes = length(unique(genes))) %>%
  # Exclude self pairs (gene in h2 and h2)
  filter(tp.x != tp.y) %>%
  # Check for duplicate pairs in reverse order (h2-h4 and h4-h2)
  mutate(comb=ifelse(tp.x < tp.y,
                     yes = paste(tp.x, tp.y),
                     no  = paste(tp.y, tp.x)))
# Remove the duplicate pairs
df <- df[which(!duplicated(df$comb)),]

And visualise it using for example a heatmap:
(ofc it looks stupid because I have a tiny sample where not all combinations are even present)
# Visualize number of genes present in each combination
library(ggplot2)
ggplot(df, aes(x=tp.x, y=tp.y, fill=N_genes)) + geom_tile()

excuse me, i did not quite understand. i have 8 separate lists of genes so how i know about 0 and 1 (absence and presence of a genes among the lists)?

@FereshTeh 0 and 1 were just placeholder values to give the info "geneY is present (1) / absent (0) at timepoint hX" . You could have expression values, ratios, pvalues or whatever variable instead of 1, and NA instead of 0, for example. Does it make sense ?

thanks a lot, actually i just have 8 lists of genes, lets say that each list has one column of genes. i want a pairwise comparison for common genes for any possible combinations of these time points.

Oh sorry, i did not understand the last part about pairwise comparisons. First, if you want to filter common genes, you could just filter rows (genes) which have values available for more than x timepoints. For pairwise comparison, maybe that'd be easier if you provided more details/background. (E.g. do you want to compute a single summary stat per condition and compute it, or is it per gene ?)

thanks a lot, i think per condition. actually i need common genes. for instance i have 630 common genes between h2 and h4 time points, while i have 1143 genes in h2 and 768 genes in h4 itself. so, i have two ways comparison; 630/1143 and 630/768. that would be great if i could do the same for all of my time points. something like having a matrix of percentage or a graph of them.

Oh I think I understand. Could you update your original question with the exact comparisons you want to make (e.g. Compute number of available genes between all possible pairs of combination). Then I can make sure I got it right and update my response.

I updated my solution, please mark it as accepted if it answers your question.

excuse me, i quite did not understand if i have 9 files each with one column of genes, how i can have df (long format). :(:(

I slightly edited my answer. hopefully it is more clear. To generate df, you need to generate the table in long format. I made a section with instructions (bash and R) to generate such table. If you chose R, you might need to edit the code to suite your data. If you chose bash, you just need to read the "long_table" file into R using read.table().

sorry for sure your code will work. i am really sorry but I'm too new in coding. for example for producing long format n r, should i have a df for each time point? i mean 9 df files? i don't know what is tp or should i put my 9 columns in  one file with its own headers for example h2, h4 in header of each column. i really could not figure out. however definitely your code is so practical

No problem with being new to programming. I am afraid I cannot help you anymore than that, but if you get a bit more into programming first, I think it will make sense. There are great free resources out there. Here is a general purpose list:https://github.com/Michael0x2a/curated-programming-resources/blob/master/resources.md
For R I would maybe try: https://www.datacamp.com/
For python in bioinformatics: http://rosalind.info/problems/locations/
For bash: https://www.codecademy.com
Hopefully some of these will be helpful to you :) Best of luck

An UpSet Plot is a useful way of visualising these kind of multi-dimensional overlaps. There's an R package for producing them, with a reasonable set of examples.
The input is a binary matrix indicating set membership:
         h2    h4    h6...
GeneA     0     1     0
GeneB     1     1     0
GeneC     0     1     1
...

And the output is a pair of bar charts, showing the size of various sets, with a dot matrix showing the set membership. An example from the UpSetR documentation:
library(UpSetR)
mutations <- read.csv(system.file("extdata", "mutations.csv",
                      package = "UpSetR"), header=T, sep = ",")
upset(mutations, sets = c("PTEN", "TP53", "EGFR", "PIK3R1", "RB1"),
      sets.bar.color = "#56B4E9", order.by = "freq",
      empty.intersections = "on")

I would simply do this...

# sepcify your timepoints, each timepoint is a list of genes!!
timepoints <- list(h2,h4,h6,h8, h10,h12,h14,h16)
# initialize an empty array that you will fill with your timepoints
vect = vector()
# go over all combinations and fill the array with the set intersection values you want (I guess...)
for (i in timepoints) {
  for (j in timepoints) {
    vect <- c(vect, length(intersect(i,j))/length(i))}}
# reshape your array to a more ordered matrix
your_matrix <- matrix(vect, nrow=length(timepoints))
# assign columns and rows names
colnames(your_matrix) <- c("h2","h4","h6","h8","h10","h12","h14","h16")
rownames(your_matrix) <- c("h2","h4","h6","h8","h10","h12","h14","h16")


text:
Identifying amino acid substitutions from local alignments in R
I want to identify the position and amino acid change in a specific region of interest in my sequences, and store that information in a table.
Is it possible to do something like this in R, by maybe using the bioconductor package? I have managed to do a simple alignment of the sequences with AlignSeqs() from the DECIPHER package, but I am unable to pull out the differences in the sequences automatically. I start with FASTA files.
I want to end up with something like this:
Isolate ID    Reference_AA   Sample_AA   Pos
1             S              T           254
2             T              D           200
3             L              A           230

I have a reference sequence which is 74 AA long, and I want to look at differences in the query sequences (which are much longer than the reference) compared to the reference, and list the results in a table. The position column is related to the position in the reference sequence and not in the query sequence. I want the first AA in the reference sequence to start at 68, not 1.
I find it difficult to add example sequences for this as they tend to be very long, but here is something shorter to work with (not related to the table above):
>reference
VGRALPDVR

>query1
KSSYLDYAMSVIVGTALPDVRDGLKPVHRRVLY

>query2
ELKSSYLDYAMSVIVGRAAPDVRDGLKPV

Expected output:

ID       Reference_AA    Sample_AA   Pos
query1   R               T           70
query2   L               A           72

What would happen if the reference is found duplicated in the queries? Should they also be returned? Do you want only 1 amino acid substitution? The last part of the first AA of the reference be 68 is easy, just add 67 to any reference position if I understood it correctly. One way would be to use regex with one position changed to ., but then one would need to work how to retrieve the position of the substitution.

If there are multiple differences in a single query, I would want each one to be listed in the output (one for each row). If there is no differences between the query and the reference, nothing needs to be written in the output

I've provided an accepted answer to this over at stackoverflow

I'm voting to close this question as off-topic because it has been cross-posted (and answered) on Stack Overflow.


text:
Aligning ChIP-Seq reads to repeats for downstream peak analysis
This is a brief question regarding the above. I have previously used bowtie to map reads from paired-end ChIP-Seq sequencing, then used the positions for peak-calling. However, I'm trying to do similar with a dataset which may consist of repeats (I have previously aligned the dataset using bowtie, and a lot of reads aligned to repeats).
For peak-calling purposes, however, I need to know reasonably confidently where precisely those repeats are located, so that MACS2 treats the coverage pileup sensibly. So my question is this: what's the best way of mapping these (short) paired-end reads from a ChIP-Seq experiment so I can do this?
EDIT: I want to be able to run a peak-caller after aligning the data, so I want a good approximation to the coverage at repeating elements after mapping for instance - I don't want to throw away all reads which multimap, as may mean that some reads could lose their mates, and so I won't be able to calculate coverage or call peaks using the data.
I want an aligner to keep both ends of the paired-end pair within about 180bp of each other (as that is the expected length of the DNA fragment that is sequenced). The immediate step after aligning the pairs is to flatten the pair into the DNA fragment that has been sequenced. With multimapped reads this step could be done by seeing which of the pair is uniquely mapped and choosing the mate which maps closest to it.
I have tried bowtie (and bowtie2) for this, which is where I understood that the dataset maps to repeats.

This is very good question, but as far as I know there really isn't a good answer. I will attempt to provide some comments and tips.

However, I'm trying to do similar with a dataset which may consist of repeats (I have previously aligned the dataset using bowtie, and a lot of reads aligned to repeats).

If this dataset is also PE and the alignment is outputting either uniquely mapped reads or a single (random) position for multimappers, it is pretty much as good as it gets. A lot of repeat regions are not exact copies, so using PE there is a good change that one read of the pair mapping to a single position is helping with mapping of the other read.

I need to know reasonably confidently where precisely those repeats are located

It is probably a matter of phrasing, but what one needs is to know the precise genomic location of the reads. The location of the repeat elements can be obtained from the Repeat Masker track of UCSC, but you probably know this already.

With multimapped reads this step could be done by seeing which of the pair is uniquely mapped and choosing the mate which maps closest to it.

I think bowtie already does this. Did you map the read pairs together or separately?
Note that:

There will always be a number of multimapping reads (fragments) that will be discarded. The information of how many are being discarded for your samples would be helpful.
Mapping can be improved using both reads and tunning bowtie settings (see bellow)

Finding a valid paired-end alignment where both mates align to
  repetitive regions of the reference can be very time-consuming. By
  default, Bowtie avoids much of this cost by imposing a limit on the
  number of "tries" it makes to match an alignment for one mate with a
  nearby alignment for the other. The default limit is 100. This causes
  bowtie to miss some valid paired-end alignments where both mates lie
  in repetitive regions, but the user may use the --pairtries or -y
  options to increase Bowtie's sensitivity as desired. source
The immediate step after aligning the pairs is to flatten the pair into the DNA fragment that has been sequenced.

As indicated above, the aligner will take care of that sort of. It will gave the mapping information of each read pair, which is then used by MACS.
Other tools
I tested ShortStack which uses the coverage of uniquely mapped reads to place multimapping reads for ChIP-seq, and the results were mixed. That said, it was SE data and the code has been through some changes, so it be worth a try.
When requested, the authors of RepEnrich provided me with the scripts for an utility called RepConsensus which maps reads to a consensus reference of repeat regions. The caveats are:

mappings are better for elements whose sequences have not diverged a lot yet
it will tell the position of the protein in the repeat sequence, but no information about the surrounding genomic regions because there isn't any.

Whatever you do, do not use all alignments, that is double counting of reads, and think very carefully before using custom methods. A blog post detailing some of the potential issues can be found here (I wrote it).


text:
Plotting metadata from GenomicRanges object
I have a GRanges object with metadata that I want to plot. Many of the online tutorials which manipulate these GRanges for plotting show how to plot the main columns (chromosome,start/end,strand) but nothing to plot the metadata as a function of the position in each row.
How precisely should I do this? Ideally, I would like to plot metadata from multiple GRanges objects on the same plot to compare them.
EDIT (following comment):
GRanges object with 6 ranges and 140 metadata columns:
      seqnames        ranges strand |
         <Rle>     <IRanges>  <Rle> |
  [1]     chr1 100000-199999      * |
  [2]     chr1 400000-499999      * |
  [3]     chr1 600000-699999      * |
  [4]     chr1 700000-799999      * |
  [5]     chr1 800000-899999      * |
  [6]     chr1 900000-999999      * |
                  eigen   domain         arm
              <numeric> <factor> <character>
  [1] -1.62475825405065     open          1p
  [2] -1.62475825405065     open          1p
  [3] -1.62475825405065     open          1p
  [4] -1.62475825405065     open          1p
  [5] -1.63354915136132     open          1p
  [6]  -1.6322733358426     open          1p
      exon.fraction       gc     short      long
          <numeric> <matrix> <numeric> <numeric>
  [1]       0.02968  0.45815         0         0
  [2]       0.04305  0.40192         0         0
  [3]       0.00959  0.42091         0         0
  [4]       0.04593  0.43337         0         1
  [5]       0.02425  0.44822         2         1
  [6]       0.16977  0.62123        27        56
         nfrags             ratio   short.corrected
      <numeric>         <numeric>         <numeric>
  [1]         0               NaN              <NA>
  [2]         0               NaN              <NA>
  [3]         0               NaN              <NA>
  [4]         1                 0 -3.50822623240367
  [5]         3                 2  2.81796238367926
  [6]        83 0.482142857142857  39.8824095079959
...

My question regards plotting a column like eigen as a function of the position (given by seqnames and ranges)

Can you show us head(., 10) for a few of the datasets so we can understand what you're trying to plot?

You can access the data.frame just using $. All of them can be accessed using values() . Something like this:
library(GenomicRanges)
gr = GRanges(seqnames=1,IRanges(start=1:10,width=1))
values(gr) = data.frame(eigen=rnorm(10),ratio=runif(10))

plot(start(gr),gr$eigen)

or
plot(start(gr),values(gr)$eigen)


text:
Tools for quality trimming at 5 prime?
I'm looking for a tool that is capable to do quality trimming at 5' end and it is configurable. For example, I can choose the quality trheshold, the read length, etc. Any recommendations?
I'm currently using sickle for quality trimming. However, one of my samples have very low quality at the beginning as well and sickle doesn't seem to be able to deal well with it.

Hi there, looking at README of sickle "You can also change the length and quality thresholds for trimming, as well as disable 5'-trimming and enable truncation of sequences with Ns." which sort of suggests that it actually trims 5'. I feel like I am missing something.

Yes, I've tried that, but I've got the feeling it's too stringent for this particular sample.

Ah, now I get it. So you search for a tool that is configurable for quality trimming of 5' end, right? I would recommend [edit] your question and explain this.

You can have a look at cutadapt. It is capable of quality trimming for both ends as you can read here.

I've just tried this, however I still see the 5' quality still very low. I run it like this:
cutadapt -q 30,30 -o output.fastq input.fastq
Is there any thing else I'm missing?

Could you please give some quality summary information of your reads? Before and after the trimming?

I never spent too much time on choosing my trimming software, therefore I might have missed some jewels. I use trimmomatic when I need versatility and I am entirely sure that it trims reads on both ends. A more convenient option is skewer, they do perform quality filtering, but I could not find an explicit confirmation that it works on 5' too.

Out of all of the major trimming tools available and widely used (trimmomatic, cutadapt/trimGalore (trimGalore is built on top of cutadapt), fastp), I actually instead prefer bbduk which is part of the bbmap suite:
https://github.com/BioInfoTools/BBMap
https://github.com/BioInfoTools/BBMap/blob/master/sh/bbduk.sh
It's much more "fine tuneable" than the other tools and runs in much less time using much less memory (as long as you set the memory limit as it's written mostly in java!). You can trim both ends or use a sliding window approach similar to what you see in Trimmomatic. Even better you can limit the number of bp in that it will trim from either left/right (5'/3'). It's by far the best all-around trimming tool for any type of sequencing dataset ({RNA, ChIP, PRO, GRO, ATAC, DNA}-seq etc...) in my opinion.


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

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

import sys
import mygene

mg = mygene.MyGeneInfo()

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

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

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

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

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

import sys
import mygene

mg = mygene.MyGeneInfo()

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

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

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

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


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

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

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

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

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

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

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

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

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

library(ggplot2)
require(gridExtra)

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

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

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

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

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


text:
Spearman correlation for large dataset
I have two datasets (DataA and DataB) and I want to find the Spearman correlation between genes and also pull out the gene names (stored in first column of dataset) in R. I am using fread from read.table to read the file and cor.test to find Rho and p-value.
The code works fine for a smaller data set of 10 obs, 10 variables. But, my actual data set is larger (DataA - 795 obj, 542 variables and DataB - 925 obs, 542 variables).
When I try to run the code on my actual data, it is consuming a large amount of time (it has been 23 hrs since my code is running). Is there a way to optimize this code? Or is there a bug in the code?
library(data.table)
A.data <- fread("DataA.csv", header = TRUE, data.table = FALSE, stringsAsFactors = FALSE)
B.data <- fread("DataB.csv", header = TRUE, data.table = FALSE, stringsAsFactors = FALSE)

correlation.result <- data.frame()

## Computing Correlation

for (i in 1:nrow(A.data)) {
  for (j in 1:nrow(B.data)) {
       correln <- cor.test(as.numeric(A.data[i,2:ncol(A.data)]), as.numeric(B.data[j,2:ncol(B.data)]), method = "spearman", exact = FALSE)
       rho <- correln$estimate
       p.val <- correln$p.value
       Gene <- paste(A.data[i,1], B.data[j,1], sep="_", collapse=" ") # To pullout gene name pairwise
       comb.data <- data.frame(rho, p.val, Gene)
       correlation.result <- rbind(correlation.result,comb.data)
 }
}

Your code looks fine, but rbind to progressively build and resize a data frame can be gluttonous in terms of memory and time, especially given the A-row x B-row times you are doing this operation, which is fine for very small datasets, but not often realistic for real-world data.
Instead, perhaps cat results to a text file, where results are written in tab-delimited columns. You can then read in that text file later on.
In addition to serializing or saving your results to text, which make them easier to work with in non-R contexts, you won't have to deal with R's memory allocation problems when running rbind repeatedly. You just write out the answer to a file and read it in later on, in one pass.
Another approach is to preallocate memory to a presized data frame, and then write to that allocated space.
You know the size of your inputs ahead of time, which makes this approach tenable. In this case, it looks like you are storing two floats and a string (rho, p.val, and Gene) and that you are going to store A-row x B-row number of them.
Here's one way you might go about this:
# create a presized data frame
Arows <- nrow(A.data)
Brows <- nrow(B.data)
correlation.result <- data.frame(t(matrix(c(0.0, 0.0, "foo"), nrow=3, ncol=Arows * Brows)), stringsAsFactors=FALSE)
colnames(correlation.result) <- c("rho", "p.val", "Gene")

# assign values to that presized data frame
Acols <- ncol(A.data)
Bcols <- ncol(B.data)
for (i in 1:Arows) {
  # get A data
  Ai <- as.numeric(A.data[i, 2:Acols])
  Aname <- A.data[i, 1]
  for (j in 1:Brows) {
    # calculate row offset in presized data frame
    ij <- (i - 1) * Brows + j
    # get B data
    Bj <- as.numeric(B.data[j, 2:Bcols])
    Bname <- B.data[j, 1]
    # run test
    correln <- cor.test(Ai, Bj, method = "spearman", exact = FALSE)
    rho <- correln$estimate
    p.val <- correln$p.value
    Gene <- paste(Aname, Bname, sep="_", collapse=" ")
    # assign value to row of result data frame
    correlation.result[ij,] <- c(rho, p.val, Gene)
  }
}

This avoids R needing to resize and recopy the data frame contents as they grow, which I suspect is the major reason it is taking so long to run.
Chapter or "Circle" 2 of Patrick Burns' book The R Inferno discusses the problems with growing objects with rbind and ways to handle it — you might find this a useful read. This book is a PDF download via your search engine of choice.

Is there a reason for the downvote(s)?

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

You can use corr.test() from psych package, e.g.
ct <- corr.test(t(A.data), t(B.data), method='spearman', adjust='none')

then ct$p and ct$r would give you (uncorrected) p-value and correlation matrices of size 795 x 925. Note that t() is required to compare matrices row-wise.
That would be orders of magnitude faster than loop-based approach.


text:
Comparing two multi-fasta files of the same set of proteins with parser - to find and count mutations after treatment
My task is to count the mutations  occurred in several proteins after a treatment. The sequences are all present in the two files in the same order. I opened both files with the FASTA parser (SeqIO.parse) in Biopython and I got all the proteins listed (separated before and after the treatment).
How can I zip the parsers together to count the mutations?
How can I count the mutations that occured after the treatment?
from Bio import SeqIO
for normal_samples in SeqIO.parse("/data/statistic/normal_samples", "fasta"):
    print(normal_samples.id)
    print(repr(normal_samples.seq))
    print(len(normal_samples))

for treated_samples in SeqIO.parse("/data/statistic/with_treatment", "fasta"):
    print(normal_samples.id)
    print(repr(normal_samples.seq))
    print(len(normal_samples))

dict_n_t = dict(zip(normal_samples & treated_samples))

from Bio import SeqIO
for normal, treated in zip(SeqIO.parse("/data/statistic/normal_samples", "fasta"),
                           SeqIO.parse("/data/statistic/with_treatment", "fasta")):
    ... do stuff...

That's generally how you zip iterators together in python.


text:
Locate KEGG/Reactome pathways in user defined network
I am working on big networks of genes. I want to locate KEGG/Reactome pathways in my network. I have already tried cytoKegg and ReactomeFI (Plugins of Cytoscape), they only give me a list of Pathways with the genes in it. But I want to visually see the grouped genes from my network which belong to one pathway in either KEGG or Reactome.
My network file looks like this
ADAM17 NOTCH1 ->
ADAM17 TNF ->
ADAM17 TRAF6 ->
ADRB2 GNAI1 ->
ADRB2 GNAI2 ->
AJUBA CDH1 <-
AJUBA CTNNB1 <-
AJUBA LATS2 -|
AMOTL2 LATS2 ->
AMPH EGF ->

I would like to see the genes in a pathway in KEGG or Reactome from my network to be either grouped together or colored. Is there any tool that does this?

How is stored your network of genes? (Or where are you working with it, you seem to be using Cytoscape but just to be sure). When do you say watch the pathways visually, do you want to get them highlighted with some color, or just focus on them, or group them aside from the main node of your network? Could you [edit] the question to add these details? Thanks

Exactly either of these, grouped or colored!

You forget to answer if you are working in Cytoscape :D. And could you add please what strategies have you looked for to find a solution? To prevent proposing the same (not-optimal) solution

I did not find any tool that could do this in R so yes I tried a few plugins from Cytoscape. But I am not confined to use cytoscape. Like I mentioned above I did try a few tools from R like KEGG itself and Pathview and plugins from Cytoscape and tried to even manually find the pathways but it gets tedious once the network starts growing.

How big is your network 1000 genes? Did you try pathview package of Bioconductor? I'm sure that in Cytoscape you can select a genes from a list and highlight them, but as I can't point to the options I'll abstain to post an answer


text:
Generating burrows-wheeler-transform and a suffix array of a DNA sequence with less memory
Also posted on biostars
After indexing
bwa index -a bwtsw reference.fa

I got files. like .bwt file and .sa file. The naive way of generating these file is:

Create all rotations of T
Sort the rotations lexicographically
Concatenate the last character of each rotation to form

BWT(T)

I used this naive way in my algorithm. The problem is it took too much memory.
My question is: is there an effective memory way to generate .bwt file and .sa file ?
Any help will be appreciated.

No, there is no way to decompress .bwt or .sa files, because they are not compressed. After indexing you'll get these files, which some programs (bwa) use. They should be in the same folder as your reference fasta file. What are you trying to do with these files? To get proper help here, please explain why you need these files, and show how you made the files.

thanks for your reply. I need to use BWT and SA in my algorithm that why I need to decompress them because in the first step to produce BWT and SA: we have to form all rotations. All these rotations took too much memory that why I prefer to produce BWT and SA by using BWA. Do you a way to use BWT and SA while they are compressed ?

Or what is the best way to produce BWT ? better than the naive version of the BWT that forms all rotations.

I still do not fully understand the problem or what you need. .bwt and .sa are needed for bwa alignment and can be generated by your bwa index command. If you run out of memory, I think the solution is to buy more RAM. Or are you talking about the size of your HD? Please ask a new question, where you exactly explain your limitations and more about your algorithm and why (and how) it needs to .bwt and .sa files.

ok here is my question in another way: how BWA generate .bwt file and .sa file in with less memory usage ? Do they save all rotations temporary ?

This is indeed a whole other question.

This is my actual question I miswrote it. Ok I will open a new question.


text:
Batch detection of CRISP proteins in fasta file
Probably a naive question. I am inexperienced.
I am interested in identifying potential CRISP (Cysteine-rich secretory proteins) in a certain tissue transcriptome (ca. 20k sequences in fasta). I have detected signalP and estimated % of cysteine in sequences. However CRISP rely on a certain pattern, and I am not aware of batch-detected algorithms.
Please would anyone know the best/simplest way of detecting CRISP candidates?
[EDIT]
I have been asked for the pattern. I summarise what the literature says:

Several of them have a CAP domain;
There should be a region called 'cysteine-rich domain' (CRD) upon the carboxyl terminal half of the protein containing 10 of roughly 16 conserved cysteines.

I am particularly interested in venom-derived CRISPs, but I am not sure whether they cluster in a certain CRISP family.
References:
Roberts et al. (2007) Structure and function of epididymal protein cysteine‐rich secretory protein‐1.
Yamazaki & Morita (2004) Structure and function of snake venom cysteine-rich secretory proteins.

Could you [edit] your question and tell us what this pattern is? Is it a specific sequence pattern? If you can define the pattern, we could be able to figure something out even if there aren't any existing dedicated tools for this.

@terdon I have read some references and hopefully the pattern described in my edit update comprises what I am looking for. Thanks

Is your goal to detect the actual proteins from the sequences? If you'll settle for detection of transcripts from the genes encoding these proteins as opposed to detection of the proteins themselves, here is an approach.
Aligning your transcriptome to a reference genome (probably with STAR) will be the first step. You have lots of short reads in your FASTA, but in their current state they are no more informative than a box of puzzle pieces. Aligning your reads to the reference will allow you to extrapolate which genes that you short reads correspond with. From here you can quantify the features (genes) in your transcriptome.
Lucky for you, there are only 3 genes in the cysteine rich secretory protein family, CRISP1, CRISP2, and CRISP3. Alignment with a reference (and annotations) and quantification of genes that your reads align to will give you all you need. You'll just have to find the three genes in your final counts table.
I recommend STAR as an aligner with the --quantMode GeneCounts argument. It's very quick if you've got the resources and is geared towards RNA-seq reads.

My fasta is the result of de novo assembly, plus we have a reference genome from a proximate species for comparison. Problem is, these are not intensively-annotated organisms (e.g. not Drosophila nor Mus musculus) thus I am trying to fish out as many annotations of novel transcripts as possible.


text:
Downloading full search results from BlastKOALA?
I have just started using BlastKOALA KEGG which has been useful in annotating (aminoacid) sequences. This is their website: https://www.kegg.jp/blastkoala/
When you get results, there are links for downloading. However these links will not download all detailed query search results, but just general notes already on the screen. To get these results I need to manually click on each query result on the page, which becomes impracticable with >500 entries. Thus I think need is a tool to download all linked contents from a webpage. I have been trying 'wget' however it doesn't work. It says 'Requested Job Not Found' whatever I do.
Please, did anyone every try to achieve this? Thanks in advance.
P.S. I have provided a working link as an example, below:
https://www.kegg.jp/kegg-bin/blastkoala_result_list?id=7345faab3969a3ef10a4a7542ba749386068138f&passwd=NKfm46&type=blastkoala&code=user
I think it should remain active for a week.
I have inserted a picture showing examples of query result links below:

At this moment I would like to end up with the listed hit results for each of the queries, whose links are showing in blue; one result per file. In short, it is using wget for downloading all links in the page into a dedicated folder.
Thanks in advance!

Could you provide some links of the entries you want to download and the command with wget you tried? This would be useful to create a program or change the command line.

@Llopis Sure, I'm adding that now to the question. I think they keep temporary links for results. Thanks!

I can't see any link to the entries you want to retry on the link you posted.  But you can download the annotations. Perhaps is that what you want?

Yes, please [edit] your question and show us an example of the data you want to download. I don't understand what you need from your question.

@terdon thanks, I have done just that. I may have found a workaround, but if anyone provides a better solution it would be great.

So you want to download the target of each of the links in the "Query" column? Do you also need the "Second best" ones? And how do you want them organized? One file per "query"? Everything in one file? Do you just want the "gene (Length)" column from the target? Please [edit] and give us an example of what you want to end up with.

Perhaps it might be easier for you to explain why your current "workaround" is not good enough for you.

@Llopis Following scattered suggestions I've gathered from related download threads I am using the command line wget -r -l3 -H -t1 -nd -N -np -erobots=off "link_path" which is returning all hits, but in a very disorganised way and raw format. I am no programmer myself.

How would it be better organized? Not being a programmer doesn't prevent you to program (or to learn to do it :) and even less to explain yourself. Please continue clarifying the question by answering terdon answers


text:
How can I improve ESS values for ucldMean in beast?
Using starBeast I ran a 50 locus tree, with 34 individuals, using evolutionary models indicated by AICc in jModel test. I ran this tree with a MCMC chain length of 500 million, storing every 10,000 generations. I used a relaxed clock log normal.
I get low ESS values for posterior, prior, birthRate.t:Species, YuleModel.t:Species, popMean, and TreeHeight.Species.
ucldMean.c:Locus_name is low for each loci.
The TreeHeight statisitcs also have low ESS values.
Also rate.c:Locus_name.mean for each loci is low as well.
Is there a way to change my parameters in Beauti to help increase the ESS values of these statistics? Is there a parameter that would likely cause the ESS values of these particular statistics to be low? Would increasing the MCMC chain length increase the ESS values, given that I've already run half a billion?

It is not my field, but I don't understand many abbreviations you use. You might get (more) answers if you explain with simpler words what is the problem. As I understand you want to increase the value of one value by optimizing some parameters of a program, is this helpful? Is there something in the mailing list of the tool helpful? Did you try or read something to solve this?

The official BEAST website suggests a number of solutions to improve the ESS of a parameter
It is not uncommon to have several billions of generations so I don't imagine increasing your number of generations would hurt.


text:
Combine multiple species trees
I have five species trees run using starBeast and Beauti. Here is a link to the paper discussing both,
http://journals.plos.org/ploscompbiol/article/file?type=supplementary&id=info:doi/10.1371/journal.pcbi.1003537.s006
 Each of these trees has the same individuals, and each locus used (10 loci per tree) has individual models set to it. All the trees used a relaxed log clock normal. I would like to be able to combine these trees into one. I have attempted to find a way to do this within Densitree and starBeast, but to no avail. When all loci are run with all trees, the program simply dies, leading me to believe that running the trees together will all parameters is computationally impossible for my hardware. Is there a way to combine these trees?

I am not familiar with trees but what are these software that you use? Could you provide a link to see what are they capable? What is a "relaxed log clock normal"? Why don't you use all 30 loci to build a single tree by combining the individual models? What do you mean that is computationally impossible: how much time do you expect to take (with how many cpus?) ?

Although I've never tried it myself, I think what you need can be achieved through the LogCombiner program from the BEAST package -- assuming the species trees files are comparable (ran with same parameters, same species names, etc.) I am also assuming that you meant 'five species tree files', each file with the posterior sample of species trees.
However, I do believe you should try to run all loci within the same analysis, or at least do some prior check to see if the distributions between data sets are comparable. LogCombiner assumes independent runs from the same data.
I would also recommend taking a look (and following) the BEAST google group (which includes starBEAST discussions).


text:
NPDock - service unavailable
NPDock service http://genesilico.pl/NPDock/ is unavailable, at least that's the message that shows after you try to submit a job. It's been like this for the last few weeks. Contact email address provided on the page also doesn't work. Does NPDock work for anyone else?

I also got the same error message; the website also says "Please run max. 5 jobs at the same time. If you need to run more jobs, please consider installing NPDock (DARS-RNP) locally.", but I couldn't find any link to the actual software.

I'm voting to close this question as off-topic because it is a question that should be asked of the tool's developers and not a general bioinformatics issue.

Sorry but we can't help with this sort of thing.

Ok, now it's fixed.


text:
Power calculation for GWAS/EWAS
I want to investigate, how much sample size i needed to obtain 80% power for GWAS/EWAS studies. Phynotype trait is discrete (not case/control) for human disease.
I wonder, does anyone has came across this? Really appreciate any suggestion or any tool/software which can answer this?

Might be worth having a look at this website: http://csg.sph.umich.edu/abecasis/gas_power_calculator/

thanks, but this if for case/control only. Our study is NOT case/control. Any other suggestion?

OK - could you edit your question to include some more details about the study you are planning (i.e. quantitative/discrete traits, species, intended outcome etc) and it might make it easier for people to give a good answer to your question. cheers.

I just Edited the question.

Do you mean continuous? If it's discrete and there are only two conditions, then the power calculations should be the same as for a case/control study.

@gringer yes I was also confused by this. what is the actual trait you are looking at?

Sorry for confusing, its a trait in young girls about their early stage in sexual maturation that happens in some higher primates and in humans, typically peaks at around 20 years of age, and is eventually involved in the development of pubic hair, body odor etc.

OK, so the trait is a continuous variable? In this case, it seems like it can be a 'time' value, like 5 years, 6 years, 5.1 years etc?

Interesting question I would recommend the effective population size $N_{e}$,  and error bars around this. Ne can be calculated via $F_{st}$ and is perhaps the most important parameter to determine the genetic population size.
Under the Fisher definition $N_{e}$ is the number of breeding individuals in the population.
The central issue with the approach you are suggesting is that population size does not equate to genetic population size. Thus the population could be huge, but if it is genetically homogeneous it will appear very small. In contrast the number of individuals could be very small, but the genetic variation between between them very large so $N_{e}$ would be huge.

Appreciated promt response.            Actually, our population  is humogeneous, and have very low sample size.  I wonder, is there any (genetics related) software where i can play around Ne and Fst?

Usual suspects Arlequin etc.. because you have low diversity the traditional stuff will work