text:
Alignment for predicting DNA hybridization?
I am currently working on a Computer Science project where we are trying to build a large set of orthogonal single-stranded DNA sequences. The goal would be to ensure that when put in solution, the ssDNA do not interact with each other to form double-stranded DNA.
I hope that this makes sense, but my knowledge  biological understanding is developing because I have a CS background. Please excuse any inadvertent errors.
As one of the first steps, we used NUPACK to predict the percent of the resulting AB complex, given two single-stranded sequences A, B in concentration 1μM each. We would use the concentration of AB to classify if the strands are likely to form double-stranded DNA.
However, we cannot use NUPACK for large sets, as we would need to do the above check against each of the other sequences, which is a quadratic problem. For a set of size 30,000 we would have to do 450 million NUPACK computations which is simply unfeasible.
So I turned to the idea of approximating NUPACK. In this paper, it is claimed that

Ideally, the thermodynamics of hybridization between a candidate
  probe and potential off-target sites would be modeled in
  silico and employed as a means of identifying probe oligos likely
  to only bind their intended targets in a given set of reaction
  conditions. Although powerful utilities such as NUPACK (52–
  54) are capable of performing such simulations, the limited
  throughput of these programs renders a direct thermodynamic
  approach impractical for genome-scale probe design. However,
  we hypothesized that features in rapidly calculated data such as
  alignment scores may be predictive of thermodynamic behavior
  and could therefore serve as a proxy for the information that
  would be produced by thermodynamic simulations.

The alignment scores are then used to predict the interaction. bowtie2 is used in the paper for alignment purposes. However, the context is different as their computation is done on genome information. Very sensitive settings of bowtie2 seem to be slow (hours) on my laptop and result in enormous files (tens of GB). I used less sensitive settings as a compromise but it seems to negatively affect the LDA.
I am looking at alternative applications for alignment, but there are many alternatives (some of the more promising ones so far, in my opinion: SNAP, MMseqs2, Genoogle). Would any of these, or others not mentioned, be better suited for my purpose?
Also any other advice or suggestions are welcome!

Welcome to the forum @David! Nice you are here and I like the question.
My question is how long (bp) are you ss probes?

Thanks for the warm welcome! We are currently interested in 20bp.

As a rough and ready and very easy to perform test I would put the entire universe of dsDNA through Primer3 (details below). This could be done by "splicing" your probe of interest onto each "chunk" of dsDNA you pass through primer3 and restrict the primer search to the 20bp primer you are interested in. In other words you are forcing the program to consider the internal and trans-hybridisation of every possible dsDNA against your probe. If you probe is not recovered from a given search you would flag this on a subsequent parsing script and then that "chunk" of DNA would be assessed further. All the components required are in the link below. You would need to provide your own wrapper, but its fairly easy.
Primer3 will assess the internal thermodynamics of your probe and it has loads of options to assess the full range of behaviour, for example 3' hybridisation (which could be extended because you want 5' to 3'). It would be easy to construct a parameter grid around all the parameters of interest. Ultimately it would provide an additional layer of certainty, and everyone knows and trusts this algorithm (its also written in C, so its fast).
The description of all the components is here , my specific purpose for customising the pipeline is different.

Bowtie2 will take a lot of RAM, I am not sure I would use it for this purpose, particularly if RAM is limiting (its used for genome assemblies against a reference). I'd ask the users whether Blat would be useful, because that is friendlier to the calculation you are performing. Anyway there some heavy duty Bowtie2 users around and their input would be useful.

Thanks for the answer! I will definitely look into it. I am a little confused though. For example, I am not sure of the meaning of "entire universe of dsDNA". Right now I am simply generating a number of random ss probes (e.g. 10,000, avoiding homopolymers) and want to make sure that any two sequences will not form dsDNA (or dsDNA is acceptable up to a threshold, say 20%). I've been testing the resulting pairs (hence the n^2) with the workflow above. In your answer I see that your are referring to a single probe, however I need to consider the entire set? Apologies if I am missing the point.

My worthwhile idea is to use Primer3 as a benchmark against your your algorithm because Primer3 is a famous tool for generating oligonucleotides. You can then assess how well your approach is against Primer3 ... we all know this program so the comparison would be meaningful to us. The basic problem with your approach is that it isn't clear how to quantify 'success' purely in silico (its a problem of most bioinformatic predictions). Anyway I think its a cool idea.

" I've been testing the resulting pairs (hence the n^2) with the workflow above. In your answer I see that your are referring to a single probe, however I need to consider the entire set? "
Primer3 will almost certainly perform that calculation, because it written in C its quick and it will select the output which best meets all the parameters. With regards your question, essentially 'yes' unless you have a heuristic search strategy. Its not hard to recursively assess all combinations and simply write increasing layers of abstraction if you have RAM problems (higher order programming)


text:
Demographic model for admixed African Americans
I am fairly new to genomics and population genetics. And I am reading this article 'Demographic history and rare allele sharing among human populations':
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3142009/pdf/pnas.201019276.pdf
In this article they propose a demographic model in Figure 4 on page 11986. However, I am confused how does African Americans fit in this demographic model that are known to have both African and European ancestry. And if I am to consider African American population than what is the split time between the two populations, time they merge into admixed population and the admixed proportion.
My apologies if it's a really basic question but, any insights would be appreciated.

Hi Shafa, welcome to the forum. Could you please provide a technical appraisal of the paper? Model "not fitting" doesn't mean very much. What quantitive aspects of the model could be wrong? "Admixture" is a good start. Generally, admixture calculations and molecular dating are quite different calculations, and admixture proportions could be intractable. It is not a basic question, if admixture meets molecular clock dating methods. Population genetic dating methods IMO are subject to inaccuracy (buts its just my opinion).

Just to further my answer, the problem with molecular dating within a species (I don't do humans [excepting immunology], but it's the same theory), is it requires an assumption of clonality (bifurcating trees). Mendelian species recombine towards panmixia, which will make any attempt to date via de novo point mutation, appear far younger than they really are, or for panmixia, incalculable (it breaks the central assumption). Dating through admixture requires a lot of extra data and/or assumptions for example the generation time,  putative admixture per generation, and shifts in population dynamics. Maybe you can calculate some of those requirements but  alot will remain assumptions.
I have not read the paper you suggested, because I'm not involved in human genetics [outside immunology], but it is possible the authors may (and I stress may) have simplified an assumption that could be difficult to estimate. It is not clear what the objective of the paper was.
Quantitative analysis is extremely important in modern population genetics/ evolutionary theory because historically it was overlooked and resulted in lots of hypothese which were difficult to test.

Thank you. I think I understand it a bit more.

source: Stack Exchange

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

sort: votes


text:
How is BLAST's nr database created?
Is there a paper or web page describing the procedure for creating the nr database used by NCBI's BLAST implementation?
I presume it's some type of clustering, but I'm curious about how exactly sequences are condensed into non-redundant representatives.

Hi, what where your search terms that didn't work? (It would help to prevent other's people time). What do you expect them to do differently than what you do?

Did you perhaps mean to leave this comment for a different question? I'm not sure what search terms you are referring to.

No, I didn't confuse your question with another. I am sure you tried to find yourself the answer to your question. I was referring to the search terms you used to find the answer to the question. I hoped we could know what did you search and suggest other search terms or explain how did we arrive to the answer (to help you in your future searches).

You are right, I did try to find the answer to my question myself. I searched for the terms in the question itself: NCBI BLAST nr, BLAST non-redundant, etc... But your comment got me to search some more and I ended up finding my answer. It's in an answer below :-)

So your answer appears to confirm my answer- good- i can delete mine if it's not helpful/needed?

Please don't delete your answer! The links to the papers are helpful and it provides more context.

Did a little more searching and found the answer in the README on BLAST's ftp site: ftp://ftp.ncbi.nlm.nih.gov/blast/db/README
6. Non-redundant defline syntax

The non-redundant databases are nr, nt and pataa. Identical sequences are
merged into one entry in these databases. To be merged two sequences must
have identical lengths and every residue at every position must be the
same.  The FASTA deflines for the different entries that belong to one
record are separated by control-A characters invisible to most
programs. In the example below both entries Q57293.1 and AAB05030.1
have the same sequence, in every respect:

>Q57293.1 RecName: Full=Fe(3+) ions import ATP-binding protein FbpC ^AAAB05030.1 afuC
[Actinobacillus pleuropneumoniae] ^AAAB17216.1 afuC [Actinobacillus pleuropneumoniae]
MNNDFLVLKNITKSFGKATVIDNLDLVIKRGTMVTLLGPSGCGKTTVLRLVAGLENPTSGQIFIDGEDVTKSSIQNRDIC
IVFQSYALFPHMSIGDNVGYGLRMQGVSNEERKQRVKEALELVDLAGFADRFVDQISGGQQQRVALARALVLKPKVLILD
EPLSNLDANLRRSMREKIRELQQRLGITSLYVTHDQTEAFAVSDEVIVMNKGTIMQKARQKIFIYDRILYSLRNFMGEST
ICDGNLNQGTVSIGDYRFPLHNAADFSVADGACLVGVRPEAIRLTATGETSQRCQIKSAVYMGNHWEIVANWNGKDVLIN
ANPDQFDPDATKAFIHFTEQGIFLLNKE

Individual sequences are now identifed simply by their accession.version.

For databases whose entries are not from official NCBI sequence databases,
such as Trace database, the gnl| convention is used. For custom databases,
this convention should be followed and the id for each sequence must be
unique, if one would like to take the advantage of indexed database,
which enables specific sequence retrieval using blastdbcmd program included
in the blast executable package.  One should refer to documents
distributed in the standalone BLAST package for more details.

Landed on that README from this question on biostars.org: https://www.biostars.org/p/217456/
Edit
In that same README file is some information on the origin of the sequences in the non-redundant sets:
+-----------------------+-----------------------------------------------------+
|File Name              | Content Description                                 |
+-----------------------+-----------------------------------------------------+
nr.gz*                  | non-redundant protein sequence database with entries
                           from GenPept, Swissprot, PIR, PDF, PDB, and RefSeq
nt.gz*                  | nucleotide sequence database, with entries from all
                          traditional divisions of GenBank, EMBL, and DDBJ;
                          excluding bulk divisions (gss, sts, pat, est, htg)
                          and wgs entries. Partially non-redundant.

Great job @juniper! It would be useful for the community if you could add to your answer that nr is a mix of different databases (not just RefSeq) and then accept your own answer. Not so humble... but that would clarify the outcome and save time for everyone I think.

Sure thing. It is done.

The Refseq team and also the NCBI resource coordinators team publish a new paper every few years, so check out the many papers (e.g. here or here), but to answer your 2nd question, non-redundancy here is (I think) defined very strictly as proteins that are identical in terms of sequence and length, so the clustering is trivial, without the need for a sophisticated clustering algorithm as required to detect more remote homologs.

This is certainly helpful reading, but RefSeq non-redundant protein sequences (a recent development, if I'm not mistaken) are not the same thing as the nr database (which has been around a while).

@DanielStandage agreed, "nr" is used ambiguously to describe non-redundant proteins (nr) or nucleotides (nt) or these Refseq records. But yes the BLAST nr databases are not just Refseq records, but also Genbank+PDB etc. i think they still use the strict identical non-redundant definition though?

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


text:
Index VCF file for rapid access by rsID?
With tabix one can index a VCF file for rapid variant retrieval based on genomic position. I'm wondering if there are any tools that will index a VCF file to allow rapid retrieval using rsIDs and/or other metadata? I'm aware of awk/grep/vcftools one-liners for this purpose, but I'd like to avoid scanning a huge VCF each time I need to retrieve the coordinates of a new batch of rsIDs.

Is re-sorting the VCF files feasible? It's possible to use Tabix on generic tab-delimited files, so sorting and indexing by the RSID column might work.

I attempted this, but it failed on most large VCFs I tried.

John Marshall said:

Moreover note that VCF files are by definition sorted by genomic position, so such an index would just be a hash table of rsID to file offset or genomic position.

I just published rsidx, a new package for indexing VCF files by rsID. As John suggests, I simply created a mapping from rsID to genomic coordinates. These are stored in an sqlite3 database.
It currently takes...several...hours to index build 151 of dbSNP on GRCh38. I'm looking into optimizing sqlite3 parameters to accelerate the indexing, as well as an alternative approach based on a two step process (populate an in-memory sqlite3 db, then dump DB to disk). But in the meantime, the search capability works great.

I am not aware of any widespread index formats for query-by-QNAME for BAM or query-by-ID for VCF. I thought there were historical samtools{-devel,-help} threads to the effect of “in principle one could index a named-sorted BAM file, but there's insufficient demand so no-one's ever implemented it”, but I can't find any just now.
Moreover note that VCF files are by definition sorted by genomic position, so such an index would just be a hash table of rsID to file offset or genomic position.
If you are using dbSNP, probably you have a dbSNP database or API that you can use to query by rsID and get genomic position back. So you can retrieve your VCF contents by rsID in two stages: do that dbSNP query, and use its results to find the variants in your VCF file by position.

high-perf-bio.
I present my high-perf-bio project, which, for the most part, solves the problem of fast extraction of anything values from VCF, BED and unformatted data.
Preparation.
1/ Run the automatic downloading of the toolkit and installation of the dependencies.
wget https://github.com/PlatonB/high-perf-bio/archive/refs/heads/master.zip && unzip -o master.zip && rm master.zip && cd high-perf-bio-master && bash __install_3rd_party.sh

2/ Reboot.
3/ cd path/to/high-perf-bio-master
Create DB with rsIDs.
python create.py -S path/to/gzipped_tables_dir -i rsIDs_field_path

Search rsIDs set in DB.
python annotate.py -S path/to/annotated_gzipped_tables_dir -D db_name -T path/to/trg_dir #-c rsIDs_column_num -f rsIDs_field_path

P.S.1. Detailed help for each high-perf-bio component: python tool_name.py -h
P.S.2. For VCF and BED you can specify the minimum of arguments.
P.S.3. Readme in English is in the early stages of writing. If you have any questions, write to Issues.

This is not too different from Daniel's project - he is more credible on the site plus has written on the  performance of his tool. Adding that information would help us compare your program to his at a glance.

@RamRS
"plus has written on the performance of his tool." The performance of such tools directly depends on the performance of the underlying DBMS. I don't think the speed difference between MongoDB and SQLite is significant for indexed data. "This is not too different from Daniel's project" My toolkit has more functionality. It contains 10 high-level command utilities for solving typical bioinformatics data retrieval tasks. Searching rsIDs is not the only feature of high-perf-bio.

Sure, your tool might be a jack of all trades but for the task at hand, is your tool faster - that is, for the same dataset, does MongoDB give you a significant advantage? In any case, my suggestion is based on the fact that you're not as well knows as Daniel, so providing at least as much info as he does will only help your case.

Take a look at grabix. Using grabix you could implement a binary search on a VCF sorted by rsID and compressed using bgzip.


text:
Database for germline copy number variations in cancer
I am interested in looking at germline copy number variations in individuals that are at high risk of developing cancer. Are there any databases where I can look, if a CNV seen in our test case has been seen/reported earlier in a patient with cancer?
Thanks.

Could you please [edit] your question and give us a few more details about what you need, exactly?

Hi @terdon, I have edited my question. Hope it clarifies the question. Thanks.

You can find CNV (and other structural variant) information in:

DECIPHER

DECIPHER is used by the clinical community to share and compare phenotypic and genotypic data. The DECIPHER database contains data from 29603 patients who have given consent for broad data-sharing; DECIPHER also supports more limited sharing via consortia. Have a look at the numbers.
Anyone can browse publicly-available patient data on DECIPHER and request to be put in contact with the responsible clinician. Why? Because sharing benefits everyone.

ExAC

The Exome Aggregation Consortium (ExAC) is a coalition of investigators seeking to aggregate and harmonize exome sequencing data from a wide variety of large-scale sequencing projects, and to make summary data available for the wider scientific community.
The data set provided on this website spans 60,706 unrelated individuals sequenced as part of various disease-specific and population genetic studies. The ExAC Principal Investigators and groups that have contributed data to the current release are listed here.

ExAC has included CNV calls since August 2016.
ClinVar

ClinVar is a freely accessible, public archive of reports of the relationships among human variations and phenotypes, with supporting evidence. ClinVar thus facilitates access to and communication about the relationships asserted between human variation and observed health status, and the history of that interpretation. ClinVar processes submissions reporting variants found in patient samples, assertions made regarding their clinical significance, information about the submitter, and other supporting data.

Database of Genomic Variants

The objective of the Database of Genomic Variants is to provide a comprehensive summary of structural variation in the human genome. We define structural variation as genomic alterations that involve segments of DNA that are larger than 50bp. The content of the database is only representing structural variation identified in healthy control samples.
The Database of Genomic Variants provides a useful catalog of control data for studies aiming to correlate genomic variation with phenotypic data. The database is continuously updated with new data from peer reviewed research studies. We always welcome suggestions and comments regarding the database from the research community.

Information from all of these is also displayed in the Structural Variant browser of VarSome*. While you can't search by patient status, so you can't use VarSome to know the clinical history of any sample with the CNV, you can at least look up a specific gene and see the structural variants listed for it. For example, for TP53, you would see:

* Disclaimer: I work for the company behind VarSome.

Thanks @terdon. I have been currently using these websites. So, good to hear this. I have also used Varsome and I think it is a nice tool.


text:
Reading a txt in R with weird header
I am trying to read a txt file in R but always I am getting
Error in scan(file = file, what = what, sep = sep, quote = quote, dec = dec,  :
  line 113 did not have 10 elements

I noticed that in the header of my file I have something like
 0000775 €   g¿´�€   4_®™00000016125 13404727020 017402  0                                                                                ustar   devons01                        stlab-icgc

Even by removing these characters I am not able to read this file
This is link of my file
https://www.dropbox.com/s/rpc5jgti5aa433b/txt.txt?dl=0
Could somebody please heap me in reading this file in R?

Where is the file from, how was it generated? Looks like an XY Problem.

Delete the corrupt portion preceeding CFS in the header and it works with with read.delim().

Sorry where is the corrupt portion preceeding CFS ?

Everything before it. Just open it in a text editor.

@DevonRyan Couldn't it be a different encoding used?

@llrs No, there was just a bunch of gibberish. It looks like some other content got added to the front of the file.

Ok, it's simply when I see this question marks I always suspect of wrong encoding...

We can skip the first row:
read.delim("txt.txt", skip = 1, header = FALSE)

Or we can get the "uncorrupted" version of the table, it looks like it is from this paper, Table S4.


text:
ComBat batch correction: understanding the model
A collaborator of mine is using ComBat for some RNA-seq data. I would like to understand what it's doing, and I have a specific question about the structure of the model. Equation 2.1 of the paper reads
$$Y_{ijg} = \alpha_g + X\beta_g + \gamma_{ig} + \delta_{ig}\epsilon_{ijg}$$
where $i$ is the sample, $j$ is the batch, $g$ is a gene, $Y$ is observed expression, $\alpha$ is an intercept term, $\beta$ is a vector of coefficients, $\gamma$ is an additive batch effect, $\delta$ is a multiplicative batch effect, and $X$ is a design matrix.
Why is $X$ not indexed by $i$? If there are $n$ samples, $m$ genes, and $d$ variables in the analysis, what are the dimensions of $X$ (or $X_i$) and $\beta_g$?


text:
Calculate the occurrence of motif in sequences (including overlaps) using Python
I need to calculate the occurrences of a motif (including overlaps) in sequences (motif is passed in the first line of standard input and FASTA sequences in subsequent lines). The sequence name starts with >, and after whitespace is just a comment about the sequence that needs to be neglected. The input of program is like:
 AT
 >seq1 Comment......
 AGGTATA
 TGGCGCC
 >seq2 Comment.....
 GGCCGGCGC
 GCCGGCGGG

The output should be:
seq1: 2
seq2: 0

I decided to save the first line as a motif, strip the comment from sequence name, join lines of sequence in one line and save sequence names (keys) and sequences (values) in a dictionary. I also wrote a function for motif_count and want to call it on dictionary values and then save it in a separate dictionary for final output. However, it does not work the way I want it as the method motif_count does not produce any result when called on my d values.
Here is my code:
 #!/usr/bin/env python3

 import sys

 sequence = sys.stdin.readlines()
 d = {}
 temp_genename = None
 temp_sequence = None
 motif = sequence[0]

def motif_count(motif, temp_sequence):
    count = 0
    start = -1
    while True:
        start = temp_sequence.find(motif, start +1)
    if start <0:
        return count
    count += 1

if sequence[1][0] != '>':

    print("ERROR")

exit(1)

for line in sequence[1:]:
    if line[0] == '>':
       temp_genename = line.split(' ')[0].strip()
       tmp_genename = temp_genename.replace('>', '')
       temp_sequence = ""

else:

      temp_sequence += line.strip()

d[tmp_genename] = temp_sequence
for key in d:
    motif_count(motif, key)

Hello Dari, this looks a bit like a homework. It's ok to ask, but it's nice to disclaim it.

Maybe it's just an error when copying your actual code in your question, but in motif_count, the while loop does not contain the test that can trigger the return, nor the incrementation of count. exit(1) should also probably be re-indented. I tried to fix this. Hopefully this corresponds to your actual code now.

Your code does nothing with the values returned by motif_count. To create a new dictionary with values modified by a function, the valmap function from the cytoolz module could be useful.

@bli please don't correct code in questions! We can't know if an error is because the OP didn't copy/paste correctly or if that's the actual code being used, so when you fix errors in the question, it could end up causing everyone to waste their time trying to find the error in code that is now correct.

I tried to fix the indentation of the code you posted, a typo in some variable names, plus a few other issues that caused it to not work as expected (see the comment in the code):
#!/usr/bin/env python3

import sys

sequence = sys.stdin.readlines()
d = {}
temp_genename = None
temp_sequence = None
# The lines returned by readlines() contain end-of-line characters,
# this is why strip() is needed
motif = sequence[0].strip()

def motif_count(motif, temp_sequence):
    count = 0
    start = -1
    while True:
        start = temp_sequence.find(motif, start +1)
        if start <0:
            return count
        count += 1

if sequence[1][0] != '>':

    print("ERROR")

    exit(1)

for line in sequence[1:]:
    if line[0] == '>':
        # Record the previous one in the dictionary
        if temp_genename is not None:
            d[temp_genename] = temp_sequence
        temp_genename = line.split(' ')[0].strip()
        temp_genename = temp_genename.replace('>', '')
        temp_sequence = ""

    else:

        temp_sequence += line.strip()

# This should be done each time a new sequence is found,
# not only at the end
d[temp_genename] = temp_sequence
for key in d:
    # We need to use the result
    # motif_count(motif, key)
    # We want to count the motif in the sequence, not in its name
    count = motif_count(motif, d[key])
    print(f"{key}: {count}")

Testing it:
$ cat tmp.fa
AT
>seq1 Comment......
AGGTATA
TGGCGCC
>seq2 Comment.....
GGCCGGCGC
GCCGGCGGG
$ cat tmp.fa | ./count_motif.py
seq1: 2
seq2: 0

Here is another way to do, using valmap to make a dict of the motif counts from the dict of sequences, and with a few other programming style changes:
#!/usr/bin/env python3

import sys
from cytoolz import valmap

motif = sys.stdin.readline().strip()

# We can make this a function of only one parameter
# because motif is fixed
def motif_count(temp_sequence):
    count = 0
    start = -1
    while True:
        start = temp_sequence.find(motif, start +1)
        if start <0:
            return count
        count += 1

seq_dict = {}
temp_genename = None
temp_sequence = None

for line in sys.stdin:
    if line[0] == '>':
        # Record the previous one in the dictionary
        if temp_genename is not None:
            seq_dict[temp_genename] = temp_sequence
        temp_genename, *_ = line[1:].strip().split()
        temp_sequence = ""
    else:
        temp_sequence += line.strip()

msg = f"""File format might be wrong"
motif: {motif}
temp_genename: {temp_genename}
temp_genename: {temp_sequence}"""

assert temp_sequence is not None and len(temp_sequence), msg

# Record the last sequence
seq_dict[temp_genename] = temp_sequence

# Make a new dict where the values are the motif counts
# valmap needs the function motif_count to take only one parameter
counts_dict = valmap(motif_count, seq_dict)
# Print the results, iterating over dict items
for (name, count) in counts_dict.items():
    print(f"{name}: {count}")

The lack of "exit (1)" indentation was an important bug

This is a good task for Biopython:
from Bio import SeqIO

with open('input_file') as f:
    motif = next(f).rstrip().lower()
    for r in SeqIO.parse(f, 'fasta'):
        print(r.id, ':', r.seq.lower().count_overlap(motif))

Note: in future I would recommend don't add the motif to the start of your FASTA file, keep it separate so as not to alter the standard file format.

The OP seems to want to read the input from stdin. Using a file is a far superior approach, of course, but maybe stdin is a requirement.

I like @Chris_Rands tidy use of OO, objects are really important per se.
There appears a paradox with your assignment that you might want to refer back to your supervisor/tutors.
If "AT" is the motif, this cannot overlap with another "AT" motif, therefore the request for "overlapping motifs" makes this part of the code superfluous.
It would be better expressed if the motif was "ATA" for example. Thus if the sequence was ATATATA the motif is present 3 times, but only twice if the motif was contiguous.
More generally the code essentially searches for AT microsatellite motifs and formally performed using e.g. tandem repeat finder.


text:
find repetitive elements on genome
I'm analyzing mygenome.fa with RepeatScout and RepeatMasker to find transposable elements. I produced a library of Repetitive elements with repeatscout and the I masked TE using Repeatmasker.
./build_lmer_table -l 14 -sequence mygenome.fa -freq mygenome.freq
./RepeatScout -sequence mygenome.fa -output mygenome_repeats.fa -freq mygenome.freq -l 14
cat mygenome_repeats.fa| ./filter-stage-1.prl >mygenome_repeats_filtered1.fasta
./RepeatMasker -s -lib mygenome_repeats_filtered1.fa mygenome.fa

Generate a masked genome using (non-low-complexity, non-tandem) repeats
cat mygenome_repeats_filtered1.fa | ./filter-stage-2.prl --cat mygenome.fa.out >mygenome_repeats_filtered2.fa

Filter out all (non-low-complexity, non-tandem) repeats that have less than 10 repeats
./RepeatMasker -pa 4 -s -lib mygenome_repeats_filtered2.fa -nolow -norna -no_is -gff mygenome.fa

After this I ran Repeatmasker:
./RepeatMasker no_is mygenome_repeats_filtered2.fa

and I produced:
mygenome_repeats_filt2.fa.masked  mygenome_repeats_filt2.fa.tbl
mygenome_repeats_filt2.fa.cat  mygenome_repeats_filt2.fa.out

I would to know a way to find the location of transposable elements masked by RepeatMasker on my genome.
Supposing that this is a part of my Repetitive elements library produced with RepeatScout as I showed above:
less mygenome_repeats_filt2.fa

>R=3 (RR=4.  TRF=0.000 NSEG=0.000)
TAAGGCGGCGAGCTGGCAGAATCGTTAGCACGCCGGGCGAAATGCTTAGCGGTATTTCGTCTGTCTTTACGTTCTGAGTT
CAAATTCCGCCGAGGTCGACTTTGCCTTTCATCCTTTCGGGGTCGATAAAATAAGTACCAGTTGAGCACTGGGGTCGATG
TAATCGACTTACCCCCTCCCCCAAAATTTCTGGCCTTGTGCCTATATTAGAAACGATTATT
>R=4 (RR=5.  TRF=0.122 NSEG=0.226)
ACACACACACACACACACACACACACATATATATATATATACATATATACGACGGGCTTCTTTCAGTTTCCGTCTACCAA
ATCCACTCACAAGGCTTTGGTCGGCCCGAGGCTATAGTAGAAGACACTTGCCCAAGGTGCCACGCAGTGGGACTGAACCC
GGAACCATGTGGTTGGTAAGCAAGCTACTTACCACACAGCCACTCCTGCGCCTATATATAT
>R=6 (RR=7.  TRF=0.134 NSEG=0.247)
TTGTTTCAGTCATTTGACTGCGGCCATGCTGGAGCACCGCCTTTAGTCGAGCAAATCGACCCCAGGACTTATTCTTTGTA
AGCCTAGTACTTATTCTATCGGTCTCTTTTGCCGAACCGCTAAGTTACGGGGACGTAAACACACCAGCATCGGTTGTCAA
GCGATGTTGGGGGGACAAACACAGACACACAAACACACACACACACATACATATATATATATATATATATA

from the file mygenome_repeats_filt2.fa.out I can see that in R=3 there is a trasposable element:
   SW   perc perc perc  query     position in query    matching        repeat              position in repeat
score   div. del. ins.  sequence  begin end   (left)   repeat          class/family      begin  end    (left)  ID

348   20.6  0.0  0.0  R=3          78   140   (81) + AmnSINE2        SINE/tRNA-Deu         67    129  (229) 234

As you can see I have the coordinates to find this element in the library, but I would to find its exactly in the file that contain my assembled genome.

Hi. In order to help you better could you explain what have you tried to do this? This way our answers can explain how to find the solution (not only the solution itself) or how to improve your skills.

ok, I modified the answer :)

Great edit!! But this code does not show what you tried to "find the location of transposable elements masked by RepeatMasker" (or at least is not clear enough for me that I am not familiar with RepeatMasker).

sorry,I'll modified again the answer and I'll try to specify it better

I modified the file with some other informations


text:
How do I add a colour annotation bar to the heatmap generated by DoHeatmap function of Seurat v2?
I am using Seurat v2 for professional reasons (I am aware of the availablity of Seurat v3). I am clustering and analysing single cell RNA seq data. How do I add a coloured annotation bar to the heatmap generated by the DoHeatmap function from Seurat v2? I want to be able to demarcate my cluster numbers on the heatmap over a coloured annotation bar. I want to have the an image like this:
 (taken from Satija V3 tutorial).
My original code:
library(Seurat)
DoHeatmap(object = object, genes.use = top20_pos$gene)

produces the following image which is not what I want:

I am aware that DoHeatmap function returns a ggplot2 object but after searching the internet for a long time, I have not been able to find  which geom layer I can add to achive my objective.
On the other hand, I have tried to generate a separate coloured annotation bar from another package called ComplexHeatmap with the intention of cropping it to fit the cluster demarcations on my heatmap generated with DoHeatmap but this has proven to be very time-consuming and inaccurrate:
library(ComplexHeatmap)

df = data.frame(type = c(rep("a", 18), rep("b", 3), rep("c", 2), rep("d", 1), rep("e", 1)))
ha = HeatmapAnnotation(df = df, col = list(type = c("a" =  "blue", "b" = "red", "c" = "green", "d" = "cyan", "e" = "magenta")))
draw(ha, 1:25)

which produces:

I will be very grateful for any hint on how to achieve my goal automatically and accurately.
Thank you in advance for your help.

@h3ab74 It is a Seurat object which has been normalised, scaled,  PCA analysed, TSNE analysed, clustered and from which differential gene expression has been calculated. All steps have been carried out with Seurat v2 methods.

bc If you can somehow derive the expression matrix being used to create this heatmap from the seurat object, you can use the ColSideColors argument in the heatmap.2 function, and pass a vector of colors through it to get this representation.. It would be something like heatmap.2(exrs.mat, ColSideColors = color_array)

If you want a simple solution, update the Seurat object to v3 and make the plot. You don't have to do the whole analysis again using v3, just use UpdateSeuratObject.

@TimStuart Thank you it worked. But is there a possibility to change the colours used on the heat map? I wish to use the blue and red colours.

You can use annotation_raster() from ggplot2 for this. https://ggplot2.tidyverse.org/reference/annotation_raster.html. If you don't need the separation (white lines in between the groups) this should be pretty straightforward. Otherwise, seurat uses empty cells (NA) to separate the groups. This is a little trickier to implement yourself.

I am not sure if it's a v3 or also v2 thing, but have you tried by setting group.by and label?
From ?Seurat::DoHeatmap:
group.by   A vector of variables to group cells by; pass 'ident' to group by cell identity classes
label      Label the cell identies above the color bar

If you prefer external tools, you can use pheatmap and setting annotation_col e annotation_colors.
First define a data.frame (df) containing the data you want in the bar (clusters in your case)
Description Condition
liver       Tumour
lung        Normal

Then you have to assign the colours
anno <- as.data.frame(df)
anno_colori <- list(
   Condition = c(Normal = '#FFC000', Tumour = '#0070C0'),
   Description = c(liver = '#D55E00', lung = '#66CC99')
)

pheatmap(mat,
     annotation_col = anno,
     annotation_colors = anno_colori,
)

The code above will add 2 bars on the top of your heatmap, one for Condition and another for Description.

There was a discussion about this in Seurat developer's github and someone actually made a package to easily add several color bars using some DoHeamap-like function. Check it here.


text:
How about BioMoose?
This is a silly post, but hope it makes a point. We have BioPerl (which never impressed me, in fact there isn't even a BioPerl tag), BioPython (which admittedly is impressive), and even now BioRust.
Why not BioMoose?
Moose is Perl5s answer to objects and is really impressive in its simplicity and ease of use. Before Moose Perl5 objects, were/are, how can I say, make evolutionary chaos theory look hip and relevant. Moose changed that completely, but stacking Moose code into BioPerl is never gonna work.

I don't understand your question. Do you ask why BioMoose is...? Or do you want to have a tag with BioMoose here on bioinformatics.SE? The reason why there isn't a BioPerl tag is that no one asked a question about it. So if you want to have the BioMoose tag, ask questions about it here.

Hi Llrs, I will provide a clearer explanation as a brief post at the end of the chat. The post is more about the future of Perl in bioinformatics, it isn't meant to be taken at face value. We don't have a tag for BioPerl because BioPerl was never a success, that was my point. If we don't have a tag about BioPython we should do (because it's good).

Sorry but I'm missing something. You will provide an explanation at the end of which chat? So, what is your question? What is the problem you are facing?

@DanielStandage summarises the status quo behind the question and is saying 'not impressed'. If you've never tried writing Perl objects (which are terrible) it's difficult emote. BioPerl likely never succeeded because it required writing Perl objects. .Moose overcomes this with aplomb, but skepticism remains.

I still don't understand what is your question (and I've briefly used BioPerl). You seem to ask why people do not use BioMoose instead of BioPerl. but that is primarily opinion-based.

Using BioPerl is not a problem, its writing it. Using objects is easy, writing Perl objects is hard and buggy. So the idea of the question is that there is a very easy way of writing Perl objects (Moose), would that enable a new Bio- repository to succeed, or (more realistically) become a Perl "revival" or "greater acceptance" within bioinformatics. I think the question is clear, and the response in one instance is "no". Perl6 might change that response and Perl5 Moose is part the transition. The foundational concepts do affect how you code.

I fear this entire thread is going to veer dangerously close to opinionation. But regardless of how one feels about the Perl language and its various libraries and communities, I think we can objectively say that the era of Perl's dominance in the bioinformatics community has passed. A lot of tools and libraries implemented in Perl are still widely used, and there is a tremendous amount of legacy Perl code around. The ability to read and write Perl remains a valuable skill.
But in 2019, few people (with a choice) are learning Perl or implementing new tools and frameworks in Perl. Python fits snuggly into the niche once occupied by Perl, and its native OO support is (I dare to say) objectively superior to any post hoc OO support Perl 3rd-party libraries provide.
The broader web/tech community moved on from Perl to Python (and Ruby and others) even before the bioinformatics community did, and this broader community's offerings (IPython, Jupyter, scipy, numpy, pytest, etc.) make Python an even more attractive alternative to Perl for bioinformatics work.
In the meantime, the "hadleyverse" ecosystem of R packages for data manipulation make R much more capable of the types of data pre-processing tasks for which an R programmer might have previously used a language like Perl or Python.
In short, Perl is waning in popularity (if not relevance), and I would not expect much enthusiasm from the bioinformatics community for a new framework built around a relatively obscure Perl extension.

Dangerously close ;-) but I definitely disagree about R, but lets see.

I will revise this question shortly and resubmit it. My perception of Python is that although the syntax is easy, precisely what the program is doing with the code appears poorly described. In Perl most coders know how the code relates to its functionality, e.g. memory allocation, where the variables reside, %INC. For example an "orchish maneuver" is understood by most coders, even if they don't know it by name.


text:
Align legends to each heatmap in a ComplexHeatmap with 2 heatmaps
I have created a ComplexHeatmap containing 2 Heatmaps, each with their own legend. The legend on the first Heatmap is continuous, whereas that of the second is discrete.
When I concatenate these Heatmaps vertically, and draw the HeatmapList, their legends are auto-aligned to the center of the plot. Is there any way I can have these legends positioned centered to the respective heatmaps and not to the center of the plot?
This is what I have:

Here's what I'd like to have:

How can I go about solving this problem? I'd really like to avoid creating customized Legends if I can, but if that's what it takes, that works too.
Thank you in advance for any help/pointers!
Note: I tried adding the ComplexHeatmap tag, but it doesn't exist yet and I don't have the reputation to create new tags yet.
Edit: Since this is a Bioconductor package, I also created an identical post on the Bioconductor support forum
Dummy data and code:
devtools::install_github('jokergoo/ComplexHeatmap');
library(ComplexHeatmap);
library(viridis);

dummy_top_mat <- structure(c(0.4563, 0.2211, 1.2235, 0.067, 1.6859, 1.0936, 0.4533,
1.6844, 1.0039, 0.5402, 0.6841, 1.498, 1.042, 0.1711, 0.3565,
2.2814, 3.0516, 1.5012, 0.5367, 3.0846, 1.1909, 0.235, 1.4266,
0.4858, 2.9054, 0.3733, 0.6902, 0.8555, 0.4234, 2.6778, 0.1568,
0.5556, 2.0172, 0.8034, 2.2897, 0.1166, 3.8033, 0.1431, 2.0606,
1.2725, 1.5365, 0.4123, 1.2087, 1.1264, 0.8334, 1.1943, 1.58,
1.5849, 0.3004, 0.3722, 0.0362, 0.0532, 1.4867, 0.4053, 0.3615,
0.0897, 1.3217, 1.1447, 1.3058, 0.1903, 0.1067, 0.9482, 1.3382,
3.2955, 0.391, 1.0418, 0.2041, 1.208, 1.5857, 3.5313, 0.472,
1.389, 0.2143, 0.0226, 0.029, 0.444, 2.0521, 0.3955, 0.3495,
0.5062, 1.3236, 1.3234, 0.7111, 0.1176, 2.2223, 1.2073, 0.3964,
2.1175, 0.3382, 0.2816, 0.71, 3.1417, 0.2402, 0.5793, 0.7662,
1.6782, 0.0986, 0.087, 0.5447, 2.6672, 1.2498, 1.0676, 1.8608,
1.8146, 0.1422, 0.4221, 0.0303, 0.9541, 0.7358, 1.7664, 1.5144,
0.2034, 0.9366, 0.7837, 0.3284, 0.1477, 1.8306, 1.3564, 0.1126,
0.0171, 2.9858, 0.0233, 0.2796, 0.6995, 1.6081, 0.215, 1.7093,
0.5178, 1.7061, 2.473, 1.8912, 0.7661, 4.4102, 1.2963, 0.6542,
0.4281, 0.4491, 0.6, 0.4076, 1.6869, 0.4747, 3.9823, 1.1226,
2.7355, 2.7036, 0.2241, 0.983, 1.0992, 1.4736, 0.1584, 0.2995,
0.2272, 0.5744, 0.9314, 0.6924, 0.0812, 1.361, 0.727, 0.1525,
1.3367, 0.566, 2.7801, 0.2349, 3.2655, 1.0675, 0.449, 0.5411,
0.8291, 0.52, 0.0507, 0.6538, 0.4636, 1.2063, 0.9784, 0.1925,
0.0756, 0.136, 1.2529, 0.317, 0.0281, 0.8668, 3.4138, 5.3898,
0.521, 0.087, 3.4819, 0.1114, 0.0061, 0.9804, 1.139, 1.4159,
1.0297, 0.6503, 1.0828, 2.3527, 0.057, 0.5602, 0.4017, 0.9985,
0.8673), .Dim = c(10L, 20L), .Dimnames = list(c("gene01", "gene02",
"gene03", "gene04", "gene05", "gene06", "gene07", "gene08", "gene09",
"gene10"), c("sample01", "sample02", "sample03", "sample04",
"sample05", "sample06", "sample07", "sample08", "sample09", "sample10",
"sample11", "sample12", "sample13", "sample14", "sample15", "sample16",
"sample17", "sample18", "sample19", "sample20")));

dummy_bottom_mat <- structure(c("Positive", "Positive", "Positive", "Positive", "Negative",
"Positive", "Positive", "Positive", "Positive", "Negative", "Positive",
"Positive", "Positive", "Positive", "Positive", "Negative", "Positive",
"Negative", "Positive", "Positive", "Positive", "Positive", "Positive",
"Positive", "Positive", "Positive", "Positive", "Positive", "Positive",
"Negative", "Positive", "Positive", "Positive", "Positive", "Positive",
"Negative", "Positive", "Positive", "Positive", "Negative", "Positive",
"Positive", "Positive", "Positive", "Positive", "Positive", "Positive",
"Positive", "Positive", "Positive", "Positive", "Positive", "Positive",
"Positive", "Positive", "Positive", "Positive", "Positive", "Positive",
"Positive", "Positive", "Positive", "Positive", "Positive", "Positive",
"Positive", "Positive", "Positive", "Positive", "Positive", "Positive",
"Positive", "Positive", "Positive", "Positive", "Positive", "Negative",
"Positive", "Positive", "Positive", "Positive", "Positive", "Positive",
"Positive", "Positive", "Positive", "Positive", "Positive", "Positive",
"Positive", "Positive", "Positive", "Positive", "Positive", "Positive",
"Positive", "Negative", "Positive", "Positive", "Positive", "Positive",
"Positive", "Positive", "Positive", "Positive", "Positive", "Positive",
"Positive", "Positive", "Positive", "Positive", "Positive", "Positive",
"Positive", "Positive", "Positive", "Positive", "Positive", "Positive",
"Positive"), .Dim = c(6L, 20L), .Dimnames = list(c("marker1",
"marker2", "marker3", "marker4", "marker5", "marker6"), c("sample01",
"sample02", "sample03", "sample04", "sample05", "sample06", "sample07",
"sample08", "sample09", "sample10", "sample11", "sample12", "sample13",
"sample14", "sample15", "sample16", "sample17", "sample18", "sample19",
"sample20")));

hmap1 <- Heatmap(dummy_top_mat, col = colorRampPalette(viridis(15))(100), name = 'log2 (TPM)', na_col = 'grey60', cluster_rows = FALSE, cluster_columns = FALSE, row_names_side = 'left', column_title_side = 'top', rect_gp = gpar(col='white', lwd=0.5), height=unit(10*0.75,'cm'), width = unit(20*0.75, 'cm'), column_title = 'Gene Expression');

hmap2 <- Heatmap(dummy_bottom_mat, col = c('Negative'='red', 'Positive'='blue'), na_col = 'grey60', row_split = factor( c('gp1','gp1','gp1','gp1','gp2','gp2'), c('gp1','gp2'), ordered = TRUE), name = 'Marker\nStatus', rect_gp = gpar(col='white',lwd=1), row_title = NULL, height = unit(6*0.75,'cm'), width = unit(20*0.75, 'cm'), row_names_side = 'left', column_title = 'Biomarkers');

draw(hmap1 %v% hmap2, ht_gap = unit(1, 'cm'), auto_adjust = FALSE);

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: macOS  10.14.4

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Users/ram/miniconda3/lib/R/lib/libRblas.dylib

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

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
 [1] viridis_0.5.1        viridisLite_0.3.0    ComplexHeatmap_2.1.0 forcats_0.4.0
 [5] stringr_1.4.0        dplyr_0.8.0.1        purrr_0.2.5          readr_1.3.1
 [9] tidyr_0.8.3          tibble_2.1.1         ggplot2_3.1.1        tidyverse_1.2.1

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.0          lubridate_1.7.4     lattice_0.20-38     circlize_0.4.6      prettyunits_1.0.2
 [6] png_0.1-7           ps_1.3.0            rprojroot_1.3-2     assertthat_0.2.1    digest_0.6.18
[11] R6_2.3.0            cellranger_1.1.0    plyr_1.8.4          backports_1.1.4     reprex_0.2.1
[16] httr_1.4.0          pillar_1.3.1        GlobalOptions_0.1.0 rlang_0.3.1         lazyeval_0.2.2
[21] readxl_1.3.1        rstudioapi_0.10     callr_3.2.0         GetoptLong_0.1.7    desc_1.2.0
[26] devtools_2.0.2      munsell_0.5.0       broom_0.5.2         compiler_3.5.1      modelr_0.1.4
[31] pkgconfig_2.0.2     pkgbuild_1.0.3      shape_1.4.4         tidyselect_0.2.5    gridExtra_2.3
[36] crayon_1.3.4        withr_2.1.2         nlme_3.1-137        jsonlite_1.6        gtable_0.3.0
[41] magrittr_1.5        scales_1.0.0        cli_1.1.0           stringi_1.4.3       remotes_2.0.4
[46] fs_1.3.1            xml2_1.2.0          generics_0.0.2      rjson_0.2.20        RColorBrewer_1.1-2
[51] tools_3.5.1         glue_1.3.0          hms_0.4.2           processx_3.3.1      pkgload_1.0.2
[56] parallel_3.5.1      clue_0.3-57         colorspace_1.4-1    cluster_2.0.7-1     sessioninfo_1.1.1
[61] rvest_0.3.3         memoise_1.1.0       haven_2.1.0         usethis_1.5.0

I think this is more suitable for GitHub issues "feature request". Currently if we add legends manually using draw(..., Legend(), ...) then we can use "just" argument within draw() to define position. But if we want to adjust existing legend within Heatmap function using Heatmap(..., heatmap_legend_param = list( it doesn't have "just" argument. (See 5.4 Section of the manual).

Even with just, it seems to be more along the lines of "where in the image" than "where with respect to the heatmap this legend belongs to". I will open a GitHub issue though, maybe the dev has a line of thought for this.

Here's the issue I opened on GitHub: https://github.com/jokergoo/ComplexHeatmap/issues/317

Turns out, this is not possible. Concatenating two heatmaps causes their legends to be grouped, and while they can be positioned together at the top, bottom, at the left or the right sides, they cannot be separated spatially.
The discussion can be viewed at the bioconductor support link from the question. Relevant snippets:

No, since you put the two heatmaps as one plot, the legends are grouped. You can either edit in Inkscape or make two heatmaps in two plots. You can get the column (order) of the first heatmap ... and assign to the second heatmap

The thing is if you add two heatmaps, the legends are always plotted together. One way I can think to improve this is to move the legends to the bottom of the heatmap.

Thanks for replying


text:
Entrez.efetch returns incomplete genbank records
I am using the biopython Entrez.efetch command to retrieve all features (CDS, mRNA, ...) of some genomes.
In this case (NC_014649, Acanthamoeba polyphaga mimivirus), it works as expected:
from Bio import Entrez, SeqIO
handle = Entrez.efetch(id='NC_014649', db='nuccore', rettype='gb')
records = SeqIO.read(handle, 'gb')
records.features

[SeqFeature(FeatureLocation(ExactPosition(0), ExactPosition(1181549), strand=1), type='source'),
 SeqFeature(FeatureLocation(ExactPosition(191), ExactPosition(2609), strand=1), type='gene'),
 SeqFeature(FeatureLocation(ExactPosition(191), ExactPosition(2609), strand=1), type='mRNA'),
 SeqFeature(FeatureLocation(ExactPosition(191), ExactPosition(2579), strand=1), type='CDS'),
...

However in this case (NC_010830, Candidatus Amoebophilus asiaticus), it only returns one record for the whole genome.
handle = Entrez.efetch(id='NC_010830', db='nuccore', rettype='gb')
records = SeqIO.read(handle, 'gb')
records.features
[SeqFeature(FeatureLocation(ExactPosition(0), ExactPosition(1884364), strand=1), type='source')]

However, if I check the corresponding file in the browser, it does have features: https://www.ncbi.nlm.nih.gov/nuccore/NC_010830.1
What could be the reason for this inconsistency ? Did I  miss something ?

This happens and when I've encountered this I resolved it by digging the features out by a different method, @arijeman 'solution' if correct needs an explanation about why because it's funndamental.

This has been answered before. See: https://stackoverflow.com/a/55402322/6262370
In short, you need to either use rettype='gbwithparts' or rettype='gb', style='withparts' to download the entire genbank flat file.

Thanks, I also had to set retmode='full', otherwise SeqIO.read would result in "ValueError: No records found in handle"

I use the efetch bash command from edirect.
With efetch -db nuccore -id CP001102.1 -format gb, I get the Candidatus Amoebophilus asiaticus 5a2, complete genome in GenBank format.
You might try your command with CP001102.1 instead of NC_010830


text:
using guppy_basecaller on node with 2 GPUs
I am trying to use the GPU enabled version of the guppy_basecaller on an HPC cluster.  I am requesting a node that has 2 GPUs and am requesting 1 of the two GPUs.  So, I am trying to set the -x cuda:-- flag.  What I was thinking this needed was the device number, ie if the allocation gave me a node and had the CUDA_VISIBLE_DEVICES= flag set, I would use -x cuda:0 or -x cuda:1 depending on which GPU I was assigned.  When I set it to 0 or auto, it seems to work when I have a whole.  When I set it to 1 when I am given a job with that device id flag set, the program errors out.  Am I thinking about the flag wrong?  The documents I have found on this are a bit sparse.
CUDA_VISIBLE_DEVICES=1
working on run id  8138
working on subset  0

git er done ...
ONT Guppy basecalling software version 3.0.3+7e7b7d0
config file:        /groups/ju_lab/softwares/ont-guppy/data/dna_r9.4.1_450bps_hac.cfg
model file:         /groups/ju_lab/softwares/ont-guppy/data/template_r9.4.1_450bps_hac.jsn
input path:         ./
save path:          ./
chunk size:         1000
chunks per runner:  1000
records per file:   4000
num basecallers:    4
gpu device:         cuda:1
kernel path:
runners per device: 4

[guppy/error] main: Could not create CUDA device 1: CUDA error at /builds/ofan/ont_core_cpp/ont_core/basecall_nn/nn_caller_cuda.
cpp:651: CUDA_ERROR_INVALID_DEVICE
[guppy/warning] main: An error occurred in the basecaller. Aborting.
```

I should have mentioned, I am trying to request 1/2 of a node because as far as I can tell guppy can't use both GPUs.  And probably should have mentioned this is on a PBS/Torque cluster using Centos 7.  And this is my command to call guppy: guppy_basecaller -i ./ --save_path ./ --flowcell FLO-MIN106 --kit SQK-LSK109 --gpu_runners_per_device 4 -r --min_qscore 7.0 -x cuda:${CUDA_VISIBLE_DEVICES} --disable_pings

New information.  Requesting a full node with 2 gpus, both -x cuda:0 and -x cuda:1 work.  But, setting the device to 1 does not result in the second gpu actually being utilized as measure by watching nvidia-smi gpu activity. So I am not sure what that flag means.  What I really need to do is specify which GPU device to use so that I don't kill other users jobs.  In SLURM with cgroups, this would not be a problem, but on this PBS/Torque cluster, I can get access to both even though I am requesting one.

OK, solved I think.  I believe -x cuda:-- is looking for the item in the list given by CUDA_VISIBLE_DEVICES rather than the actual device id.

Welcome to the site. You can [edit] the question to include the relevant code into the question. Also to mark as solved you can (and are encouraged) to post your solution as an answer (and then mark it as the accepted one). This way others, like me we, can learn from you :D

Did using -x auto not find the device you wanted?

please where did you get the software of guppy_basecaller fro Centos 7? thanks


text:
How to learn/get started bioinformatics?
I know this is a broad question and I'm sorry if this is not the right place to post (this is my first post ever), but is there some sort of roadmap for bioinformatics?
I'm an undergraduate student and I have a strong background in statistics, computer science, and molecular biology so I feel like I should be able to get around, but every time I talk or read about bioinformatics I feel lost. It's such a vast area that I don't even know where to begin or what even exists.
Is there some kind of all-encompassing textbook or roadmap for me to follow?
Most textbooks I could find expend way to much time explaining really basic stuff like what is a cell and how to program and I feel they fail to convey the bigger picture of bioinformatics. Things like what exists in the area, what are the current technologies, conventions, and challenges that we face. What should I know in order to call myself a bioinformatician?

Welcome to the forum @B-Aral. How good is your Python? A good programmer (eg. could you implement method chaining?) who is keen on extending their biological understanding is a fantastic start. However, if you are statistics orientated R is a better option. Thereafter you need to pick an area and hone your skills within that area, Python vs. R will shape which area on some instances. It doesn't matter if you're projects, post-grad training don't fit, core R or Python skills are the primary skill set IMO.

Indeed, This is an extremely broad question. There would be great disagreement about what even bioinformatics means among people in this community and I think that the only common denominator of all of us is at least a basic knowledge of programming. So, is there a way how narrow your question? Why do you want to learn bioinformatics? Are you interested in sequence analysis (the largest bioinfo community, but certainly not the only)?

Have you looked at other sites? Like this post with links to other relevant discussions but you can find some other posts by searching around.

I find Rosalind (http://rosalind.info/problems/locations/) a great resource, which in addition to programming exercises will explain biology and common bioinformatics algorithms.

The formal answer is O'Reilly "Bioinformatics programming, using Python". Check out the classic O'Reilly texts on "Mastering Perl for bioinformatics", and possibly"Beginning Perl for bioinformatics" (which is likely far too basic), which establish some core principles and have a solid reputation in teaching. However, the Perl/Python thing might get confusing (I dunno).
If you are statistics orientated at a guess, O'Reilly's, Bioinformatics data skills... I don't know, but something that would focus on R.
The real goal is core-technical competance because the biological understanding will be taught by your supervisor. Thereafter you might turn your attention to teaming up with a lab for project training.

I'd advise against learning Perl for bio-informatics these days. The only places I've run into Perl in bio-informatics is in legacy code.

@Pallie you really shouldn't ever advise anyone against learning something. Especially not something as useful as Perl. Learning is always good. In any case, while Perl is indeed not as popular as it used to be, it's still around (have a look at the Ensemble API, for instance), and it's still useful. And the specific language used is only a minor detail: whatever language you learn, you need to learn the same general concepts.


text:
What do negative dN or dS mean in cal_dn_ds in Biopython?
I'm noobie in bioinformatics. I use cal_dn_ds from Bio.codonalign.codonseq with method 'NG86'. Sometimes I get dS=-1. What does it mean? How can number of substitutions be negative? Does it mean that something is wrong with my data?
Upd.
1) Chris_Rands pointed to sources, I read them already, but the question remain. What is that ps variable, and why if it's greater than 3/4, dS = -1?
2) Some more details.
cal_dn_ds requires sequences from CodonAlignment object, which is created by function Bio.codonalign.build that requires (quote from docs) protein alignment and corresponding nucleotides. I have nucleotide sequences from NCBI at the beginning. So I trim non-coding regions from nucleotide sequences, translate all trimmed sequences, write them to a file, perform multiple alignment of the protein sequences, and then pass the list of aligned protein sequences and the list of nucleotide sequences to the build function.
I analize protein DKK1, and example sequences that get negatives are Empidonax traillii XM_027898752.1 and Physeter catodon XM_007110600.1. But I just want to understand general idea behind dS = -1, what can it mean it different cases.
3) Michael G. gave an interesting suggestion that it means zero, but it's strange because the function returns dN and dS separately, so Biopython might as well return zero, so the user would handle those cases himself. And for some sequences Biopython actually returns normal positive dS and dN = 0 at the same time. If it's an error, what are the possible causes? And why in those error cases Biopython can calculate dN, but dS is -1?

Please [edit] your question and give us some more details. What kind of input do you have? What are you passing as input to cal_dn_ds? If you are using some sort of alignment, what species are involved? The more details you give us, the better your chances of getting a useful answer.

https://github.com/biopython/biopython/blob/master/Bio/codonalign/codonseq.py#L390

terdon, Chirs_Rands, thank you for replies, I updated my question with more details.

dS = -1
.. there are two suggestions

I suggest its an error, possibly based on nucleotide calls within the S, check for anomylous nucleotides, particularly over a non-synonymous codon
It might mean their are no synonymous mutations. It is not the best idea to assign this a value of 0 in a ratio because it will cause the ratio to fall over when in fact their could be strong selection pressure, ie non-synon mutations. Having an incrementally small value is problematic because dN if positive could anything up to infinity, so -1 might flag this problem.

Okay I now know why... You are correct to be concerned with the result because its serious.
My sincere advice to you is do not use any dN/dS ratio where -1 is occurring under any circumstances. On this occasion I'm not going to explain why for two reasons:

First and foremost, as a phylogeneticist you should know why, it relates to the fundamentals of the discipline.
(Hint) The answer is in my past posts on this forum.

To conclude, I'm pretty nuanced in statistical dogma, you can check my record, IMO most dogma is tribalism because every calculation has weaknesses. HERE HOWEVER, there is a real problem and again I urge you to ditch any data set where -1 is occurring and everything associated with it. That is a fairly dire warning should give you a clue about whats wrong. If you are not prepared to do this, then you must report the result and the situations where it is arising (don't hide it).

Thank you for your reply, I provided some more details for you at Upd 3. I checked my nucleotide sequences, and they do not have anything except ATGC, if that's what you meant by 1.


text:
Multi-line efetch digging
I'm using efetch to dig out a middle "envelope" protein within a viral genome specifically within fasta format. The code below works but may not be robust for an entire database. Thus I want to replace the two "perl" lines with formal code within the Entrez package. I tried using efilter -query but couldn't get it to work. Genome positions are one approach but are not absolute from one virus to another. Rewriting this may be the best option.
ncbi2=NC_002021
efetch -db protein -format gpc -id "$ncbi2" |\
xtract -pattern INSDSeq -element INSDQualifier_value |\
perl -p -e 's/^(.+?)\s+taxon.*.*\s+\b([A-Z]{300}[A-Z]+).*/>$1\n$2/'|\
perl -p -e "s/[\t \;\/]+/_/g"

Output
>InfluenzaAvirus_A_PuertoRico_8_1934
MDVNPTLLFLKVPAQNAISTTFPYTGD.....

I will provide a more detailed explanation shortly.
Any ideas welcome.

In the meantime, you could streamline the perl a little to | perl -pe 's/^(.+?)\s+taxon.*envelope protein.*\s+\b([A-Z]{300}[A-Z]+).*nonstructural protein 1.*/>$1\n$2/; s/[\t ;\/]+/_/g. There's no need to process the file twice and you also don't need to escape the ; inside a character class.

Hi @terdon, yes and thanks now you've mentioned it, it's obvious.

Not exactly what you were looking for but may be you can use something like this? efetch -db protein -format gpc -id "$ncbi2" | xtract -pattern INSDSeq_feature-table -group INSDFeature -if INSDFeature_key -equals 'mat_peptide' -tab '\n' -element INSDFeature_location,INSDQualifier_value | grep 'envelope'

Perfect it works a treat @vkkodali. Huge thanks. A small amount of regex will finish the fasta format. Maybe you would like to provide a formal answer so I can sign this question as complete?

Not exactly what you were looking for but may be you can use something like this?
ncbi2=NC_002021
efetch -db protein -format gpc -id "$ncbi2" \
    | xtract -pattern INSDSeq_feature-table \
        -group INSDFeature \
        -if INSDFeature_key -equals 'mat_peptide' \
        -tab '\n' \
        -element INSDFeature_location,INSDQualifier_value \
    | grep 'Segment 2'

For explorations of the INSDSet XMLs, xtract has -insd option that you can use. In this particular case, you can do something like this:
efetch -db nuccore -id "$ncbi2" -format gpc \
    | xtract -insd mat_peptide product peptide

For more information, checkout EDirect documentation and some  example queries

Thank you its perfect. Stacked together in a small subroutine and parsed with a bit of regex its fantastic. "esearch -query" is used to obtain the accession numbers.


text:
Remove variants that do not map to human genome
[This question was also asked on Biostars]
I received an hg38 VCF file that's had variants imputed with 1000 genomes. I've encountered some issues with the VCF; REF alleles that do not align to a reference genome, ALT alleles that do not appear to be reported anywhere in the literature, and, most recently, variants that flat-out do not align to the human genome (variants on chr19 with bp-pos 100 million+ when the whole chromosome is in the 50 million bp range).
I've worked out hack-y solutions to most of the issues that I've encountered, but this latest one has been an issue for me. I only detected these variants when I ran VEP and it flagged them as not mapping to the genome. As such, I'm more or less removing these variants one at a time using grep -v. I'd like a solution where I can just remove any variants from the vcf that appear to map to regions that do not exist in the human genome. Bonus points if the solution also encompasses some of the other issues I mentioned, although I think I've already found solutions to those. Is there anything out there that does this?

"ALT alleles that do not appear to be reported anywhere in the literature": this isn't necessarily a problem, of course. It's entirely likely that you will find novel variants in any real life sample.

While I do agree with you, I've seen enough issues with the VCF that I'm skeptical of anything 'novel' reported in there. I've taken the conservative route and removed these variants, at least until I can get the original genotyped calls from our collaborator and verify that these calls were truly there.

Oh sure, depending on what data you expect it might be an issue. I only pointed out it doesn't have to be one. For example, I work with clinical "real world" samples and I expect several hundred, or thousand novel variants with no frequency data, let alone any mention in the literature.

I would return this VCF to whoever gave it to you and ask them why the positions on chromosome 19 aren’t anywhere near the reference positions. There is virtually no point in proceeding with any kind of analysis if you don’t even know whether the positions are correct.

Answer I received on Biostars that appears to work as expected:
bcftools norm -Ob --check-ref wx -f <reference_fasta> <vcf.gz> -o <output_name>

I would throw away this VCF. The command line only fixes obvious issues but you can never be sure to fix everything without understanding the error source. Those remaining issues may bite you and lead to false conclusions.

While in general I would agree with you, I need the vcf to get started on our analysis. We're working on getting the genotype array calls vcf that does not have the imputed data, and hopefully these issues will not exist in that vcf.

That chr19 example indicates something deeply wrong with the VCF. This could be due to inconsistent genome versions or low-level bioinformatics errors. The frequency, the alleles and even the positions could be all messed up. bcftools can't fix these for you. You may need to redo all analyses again. Ignoring the errors is not saving time; it is wasting time. Anyway, it is your call.


text:
bowtie2 options when mapping stranded single end reads
I'm trying to map RNA-Seq reads generated using the NEB Ultra Directional kit is first strand reversed protocal. I'm using bowtie-2.3.5_1 to align against the mus musculuous (mm10) genome.  The reads are single ended.
It seems that --fr and --rf are relevant only for paired end options. I tried these options as well as not specifying either (default).
E.g.
$ bowtie2 --fr \
        --threads 20 \
        --very-sensitive \
        -x /reference/mus_musculus/mm10/ucsc_assmebly/illumina_download/Sequence/Bowtie2Index/genome \
        -U trimmed_out/my_trimmed.fq.gz \
        -S bowtie_fr_1.sam

$ samtools sort -n -o bowtie_fr_sort_1.sam bowtie_def_1.sam

I then compared the output of bowtie_rf_sort.sam bowtie_default_sort.sam and bowtie_fr_sort.sam.  They are basically identical excluding the header.
Question:
Does strandedness matter in this alignment? Perhaps I should care about the --nofw/--norc options? Am I missing something here?

Are you taking transcripts as the reference? --norc might be what you are looking for.

I'm using mm10 as my reference (downloaded from Illumina).

Then you need an RNA-seq mapper, not bowtie2. STAR and HISAT2 should have options to enable strand specific mapping. I am not familiar with those, though.

Could you please elaborate how bowtie is not an rna-seq mapper? Per their website "Bowtie 2 is often the first step in pipelines for comparative genomics, including for variation calling, ChIP-seq, RNA-seq,"

Bowtie2 is not aware of splicing. For typical RNA-seq applications, you will want to use a splice-aware mapper, such as star and hisat2, which is specifically designed for RNA-seq. Tophat2 does use bowtie2 for mapping, but it is invoking bowtie2 in a non-standard way and is generally thought to be superseded by star and hisat2 anyway.

I repeat what user172818 says...don't waste your time on bowtie or tophat.  Use STAR or HISAT2, or a pseudoaligner like kallisto with a transcriptome.


text:
Determine reference for reference-compressed SRA file
Note: this question was also asked on Github.
I have 241 SRA files that appear to be reference compressed. I didn't even know this was a thing until I tried to convert them to Fastq files without an internet connection. I got the "name not found while resolving..." error described at the bottom of this page when I attempted to do the conversion.
The instructions say:

If you have elected to prevent the toolkit from contacting NCBI, you will need to manually acquire the reference(s) here

We have many NCBI databases mirrored locally, including RefSeq, so it's likely that we already have all the references. I have two questions:

How do I determine the reference used to compress a reference-compressed SRA file? (Note this needs to be scripted for 241 SRA files, so manual inspection is out of the question.)
Once I've determined the correct reference and ensured it is available locally, how to I tell fastq-dump where to find it?

Run align-info  to see the references used.
Run prefetch </path/to/sra/file> to download missed references.


text:
blast command line
I would to use blastx from command line to analize a library produced with RepeatScout in format .fasta to find new transposable elements.
I use this command mylaptop$ tblastx OB150DEC_2_repeats.fa
but it reports an error:
USAGE
  tblastx [-h] [-help] [-import_search_strategy filename]
    [-export_search_strategy filename] [-db database_name]
    [-dbsize num_letters] [-gilist filename] [-seqidlist filename]
    [-negative_gilist filename] [-entrez_query entrez_query]
    [-db_soft_mask filtering_algorithm] [-db_hard_mask filtering_algorithm]
    [-subject subject_input_file] [-subject_loc range] [-query input_file]
    [-out output_file] [-evalue evalue] [-word_size int_value]
    [-qcov_hsp_perc float_value] [-max_hsps int_value]
    [-xdrop_ungap float_value] [-searchsp int_value] [-sum_stats bool_value]
    [-max_intron_length length] [-seg SEG_options]
    [-soft_masking soft_masking] [-matrix matrix_name]
    [-threshold float_value] [-culling_limit int_value]
    [-best_hit_overhang float_value] [-best_hit_score_edge float_value]
    [-window_size int_value] [-lcase_masking] [-query_loc range]
    [-strand strand] [-parse_deflines] [-query_gencode int_value]
    [-db_gencode int_value] [-outfmt format] [-show_gis]
    [-num_descriptions int_value] [-num_alignments int_value]
    [-line_length line_length] [-html] [-max_target_seqs num_sequences]
    [-num_threads int_value] [-remote] [-version]

DESCRIPTION
   Translated Query-Translated Subject BLAST 2.5.0+

Use '-help' to print detailed descriptions of command line arguments
========================================================================

Error: Too many positional arguments (1), the offending value: OB150DEC_2_repeats.fa
Error:  (CArgException::eSynopsis) Too many positional arguments (1), the offending value: OB150DEC_2_repeats.fa

I don't understand what does it mean too many positional arguments

A positional argument is one that doesn't follow something like -foo or --bar. That is, the variable that it should internally map to depends on it's position in the command. You're looking for something like tblastx -db something -query OB150DEC_2_repeats.fa

ok. I have to download a database to use?

Or point it to a remote one.

I've never used blast from command line before. Could you explain to me how to point it from remote one? I'm working with a cephalopoda genome

Just use nt and I think it will figure it out. We only use local databases, so if that doesn't work you'll need to read the documentation.

ok. I've tried but the message is BLAST Database error: No alias or index file found for nucleotide database [nt] in search path [/mydirectory/OB150DEC::]

Use the -remote option https://www.ncbi.nlm.nih.gov/books/NBK279680/

Thank you for your help, but it doesn't work. It reminds a lot of errors. I'll try downloading an nt database


text:
Which values from Seurat::FetchData function are to be used for correlation analysis between genes?
I want to perform a correlation test between genes in on my single cell RNA seq data set. I perfomed the differential expression analysis using the Seurat version 2 package, after performing stages of normalisation, scaling, PCA, TSNE analyses and clustering. My correlation plot is far from what I expect and I want to know if the expression values obtained with Seurats V2's FetchData function (imputed, scaled, or raw) may have an effect on correlation analysis' outcome. This is what I did:
# load Seurat v2
library(Seurat)

# load corrplot package
library(corrplot)

# obtain expression values of the genes I want to test rom the Seurat object
expression_matrix <- FetchData(object, vars.all = c("gene1", "gene2", "gene3"))

expression_matrix <- round(expression_matrix , 1)

cor_expression_matrix <- cor(expression_matrix)

corrplot(cor_expression_matrix , is.corr = T, method = "square")

Which gene expression values should be best used for the analysis of correlation between genes: imputed, scaled, or raw?
Thank you in advance for your kind response.

You should use normalized data. The function has a parameter slot="data" that let's you select to access the normalized data.
Any reason for rounding to 1 decimal place here? I imagine that would have a substantial impact on the correlation you compute.


text:
bcl2fastq2 v2.20.0 converstion software on Ubuntu 18.04
bcl2fastq is a program that converts Illumina run output from .bcl to fastq. I'm having difficulties installing the software I downloaded from Illumina on my Ubuntu 18.04. I've followed their installation manual, but the software is not building when I enter the command "make".  I receive the error message... "make: *** No targets specified and no makefile found. Stop. Has anyone had similar experience? Is there another software that will install easier on Ubuntu?
Software:
https://support.illumina.com/downloads/bcl2fastq-conversion-software-v2-20.html
Manual:
https://support.illumina.com/content/dam/illumina-support/documents/documentation/software_documentation/bcl2fastq/bcl2fastq2-v2-20-software-guide-15051736-03.pdf
Steps followed:
export TMP=/tmp
export SOURCE=${TMP}/bcl2fastq
export BUILD=${TMP}/bcl2fastq2-v2.20.0-build
export INSTALL_DIR=/usr/local/bcl2fastq2-v2.20.0

cd ${TMP}
tar -xvzf bcl2fastq2-v2.20.0.499-Source.tar.gz

mkdir ${BUILD}
cd ${BUILD}
chmod ugo+x ${SOURCE}/src/configure
chmod ugo+x ${SOURCE}/src/cmake/bootstrap/installCmake.sh
${SOURCE}/src/configure --prefix=${INSTALL_DIR}

cd ${BUILD}
make

The final command is supposed to be make install, but the configure step fails with:
[ . . . ]
-- Configuring incomplete, errors occurred!
Couldn't configure the project:
/tmp/bcl2fastq2-v2.20.0-build/bootstrap_cmake/bin/cmake -H'/tmp/bcl2fastq/src' -B'/tmp/bcl2fastq2-v2.20.0-build' -G'Unix Makefiles' -DBCL2FASTQ_NAME_SHORT:STRING="bcl2fastq" -DBCL2FASTQ_NAME_LONG:STRING="BCL to FASTQ file converter" -DBCL2FASTQ_COPYRIGHT:STRING="Copyright (c) 2007-2017 Illumina, Inc." -DBCL2FASTQ_VERSION_MAJOR:STRING="2" -DBCL2FASTQ_VERSION_MINOR:STRING="20" -DBCL2FASTQ_VERSION_PATCH:STRING="0" -DBCL2FASTQ_VERSION_BUILD:STRING="422" -DBCL2FASTQ_VERSION:STRING="2.20.0.422" -DBCL2FASTQ_SOURCE_DIR:STRING="/tmp/bcl2fastq/src" -DBCL2FASTQ_PREFIX:PATH="/usr/local/bcl2fastq2-v2.20.0" -DBCL2FASTQ_EXEC_PREFIX:PATH="" -DCMAKE_INSTALL_PREFIX:PATH="/usr/local/bcl2fastq2-v2.20.0" -DBCL2FASTQ_BINDIR:PATH="" -DBCL2FASTQ_LIBDIR:PATH="" -DBCL2FASTQ_LIBEXECDIR:PATH="" -DBCL2FASTQ_INCLUDEDIR:PATH="" -DBCL2FASTQ_DATADIR:PATH="" -DBCL2FASTQ_DOCDIR:PATH="" -DBCL2FASTQ_MANDIR:PATH="" -DCMAKE_BUILD_TYPE:STRING="RelWithDebInfo" -DCMAKE_PARALLEL:STRING="1"
Moving CMakeCache.txt to CMakeCache.txt.removed

Please [edit] your question and add: i) the URL of the software you downloaded; ii) the exactly commands you have run. Are you running make from the directory where you have run .configure? Did you run configure?

Thanks for the edit, that's clearer! Now, some of those commands should have given an error message, did they? Especially the ${SOURCE}/src/configure --prefix=${INSTALL_DIR}, did that finish successfully? And what happens if you now run cd /tmp/bcl2fastq2-v2.20.0-build; make? Have you installed the build-essentials (sudo apt install build-essentials) package on Ubuntu?

No errors displayed when I ran the commands.  I haven't been able to install the build-essentials... I get the message "Unable to locate package build-essentials.  The message I get when I run "make" is something about not having a target.

please run configure again and show us the output. There should be a message there telling you what's missing. And sorry, that should have been sudo apt  install build-essential, no final s.

I'm seeing the errors at the bottom.  I'm attaching a screenshot.  The screenshot is too long, so I'll add them in three parts.

Still too long, so I'll attach a txt file.

Let us continue this discussion in chat.

Just add the last few lines, that's all we need. And did you install build essential?

I did add the build-essential

So, I started the process over and I get the error when I run "make"
laparic1@crb-byphhv2:/tmp/bcl2fastq2-v2.20.0-build$ cd ${BUILD}
laparic1@crb-byphhv2:/tmp/bcl2fastq2-v2.20.0-build$ make
make: *** No targets specified and no makefile found.  Stop.
I wonder if I the command before it, cd ${BUILD}, is correct because I'm already in the BUILD directory.

Please [edit] your question and add all this information. You can use the formatting tools to format your commands and their output as code. Comments are hard to read, easy to miss and can be deleted without warning.

I suggest you install it via the RPM, this is an incredibly annoying package to compile yourself (this is unfortunately common with software from Illumina).

Thanks. It certainly looks that way.

Install it with conda.
conda install -c dranew bcl2fastq

Hope you don't get sued by Illumina for illegal distribution of their software...

... need to preferably install Anaconda3/minicoda .

Where are you running the software? In general and especially for servers and high performance clusters I would recommend to use a module system like Lmod in conjunction with EasyBuild, which contains configs for bcl2fastq2 that take care of the installation. A module system allows you to use different versions and to keep your environment clean among other benefits.

The OP is running the software on their own Ubuntu machine, I'm not sure EasyBuild would be appropriate in context.


text:
RNASeq read coverage in protein space?

I used samtools depth to find the per base-pair read coverage over a number of isoform contigs from my Trinity assembly.
I have also conducted a multiple sequence alignment of those isoforms using Clustal Omega DNA aligner.

In my final step, I want to assess the per-residue read coverage of the consensus amino acid sequence to demonstrate an idea of confidence in that sequence.
Does anyone have any tips on how to go about doing this and whether there are existing tools for such things?

Answer from @cephbirk, converted from comment:
What I ended up doing was to choose the read coverage for the lowest covered nucleotide in each codon as each residue's read coverage.


text:
Detecting broad peaks in sRNA-seq data
What kind of tool would be appropriate do detect "broad peaks" in small RNA-seq sequencing data?
MACS2 appears to be  developed for ChIP-seq data, but I see that there is a --nomodel option. Would that make this program usable in my case?
I suppose I should also set the --shift and --extsize options, but I'm not sure how:
  --nomodel             Whether or not to build the shifting model. If True,
                        MACS will not build model. by default it means
                        shifting size = 100, try to set extsize to change it.
                        DEFAULT: False
  --shift SHIFT         (NOT the legacy --shiftsize option!) The arbitrary
                        shift in bp. Use discretion while setting it other
                        than default value. When NOMODEL is set, MACS will use
                        this value to move cutting ends (5') towards 5'->3'
                        direction then apply EXTSIZE to extend them to
                        fragments. When this value is negative, ends will be
                        moved toward 3'->5' direction. Recommended to keep it
                        as default 0 for ChIP-Seq datasets, or -1 * half of
                        EXTSIZE together with EXTSIZE option for detecting
                        enriched cutting loci such as certain DNAseI-Seq
                        datasets. Note, you can't set values other than 0 if
                        format is BAMPE for paired-end data. DEFAULT: 0.
  --extsize EXTSIZE     The arbitrary extension size in bp. When nomodel is
                        true, MACS will use this value as fragment size to
                        extend each read towards 3' end, then pile them up.
                        It's exactly twice the number of obsolete SHIFTSIZE.
                        In previous language, each read is moved 5'->3'
                        direction to middle of fragment by 1/2 d, then
                        extended to both direction with 1/2 d. This is
                        equivalent to say each read is extended towards 5'->3'
                        into a d size fragment. DEFAULT: 200. EXTSIZE and
                        SHIFT can be combined when necessary. Check SHIFT
                        option

In my case the tags are the same things as the fragments, so this "extension" thing is probably not appropriate.
Basically, I would like to call peaks based on the plain real coverage of my small reads. Are there other tools that would be more appropriate?

Do you actually need peaks or is your end goal instead to find novel sRNAs? I would think stringTie followed by a considerable amount of filtering would be better.

I mapped sRNA-seq reads on the genome, and I want to automatically delineate genomic regions that produce them. By "broad peaks", I mean such regions. Not sure what you mean by "novel sRNAs".

sRNAs shouldn’t appear as broad peaks unless they’re coming from repeat clusters. See the answer from Sebastian Müller for good tool options.

The sure why any software should be confined to only broad peaks, but there is two packages that should do what I thing you have in mind. Just to clarify, for smallRNAs peak detection is sometimes also referred to as loci or cluster detection:

ShortStack which is developed for alignment, annotation, and quantification of small RNAs. The de-novo cluster finding option as copied from the manual might do the trick:

DE-NOVO CLUSTER FINDING
    Unless options --locus or --locifile are used (see below), ShortStack
    will de-novo identify clusters of small RNA accumulation genome-wide.
    Cluster definition is simple: First, all regions containing at least one
    primary alignment are found where the maximum distance between the ends
    of the alignments is <= option --pad (default: 75). Second, if the
    number of alignments in the cluster is >= option --mincov (default:
    0.5rpm), the cluster is kept. The mincov threshold can also be specified
    in terms of reads per million by using a value such as 3.2rpm (which
    specifies the threshold to be 3.2 reads per million). Using a rpm
    threshold allows the sensitivity of cluster discovery to be normalized
    between libraries of different sizes. Alternatively, reads per million
    mapped (rpmm) can be specified: A --mincov of 1.2rpmm indicates 1.2
    reads per million mapped is the threshold. rpm is a fraction of total
    library size, while rpmm is a fraction of only the aligned & placed
    fraction of the library.

segmentSeq
is an R-package useing a Bayesian approach for small RNA locus detection. It's however really slow and not developed any further by the main author.


text:
Error: sort order of reads in BAMs must be the same
I am a novice to bioinformatics trying to run a cuffdiff job in cufflinks and keep getting this error message:
Error: sort order of reads in BAMs must be the same

GENERATING SAM FILES:
cd /fs/project/PMIU0157/Claire_Stoll/RNASeq/B73_hisat2/B73_hisat2_11-5
module load hisat2
hisat2 -k 1 -p 8 -S 11-5_B73_accepted_hits.sam -x /fs/project/PMIU0157/Claire_Stoll/RNASeq/B73_hisat2/B73_Zea_mays.AGPv4.dna.toplevel -1 /fs/project/PMIU0157/Claire_Stoll/RNASeq/11-5_FCHHL7JBBXX_L7_DLPRMzP300H3AAAPEI-27_1.fq  -2 /fs/project/PMIU0157/Claire_Stoll/RNASeq/11-5_FCHHL7JBBXX_L7_DLPRMzP300H3AAAPEI-27_2.fq

cd /fs/project/PMIU0157/Claire_Stoll/RNASeq/Mo17_hisat2/Mo17_hisat2_11R-5
module load hisat2
hisat2 -k 1 -p 8 -S 11-5_Mo17_accepted_hits.sam -x /fs/project/PMIU0157/Claire_Stoll/RNASeq/Mo17_hisat2//Mol17_Ref_Lai -1 /fs/project/PMIU0157/Claire_Stoll/RNASeq/11R-5_DLPRMzP300B4AAAPEI-2_1.fq  -2 /fs/project/PMIU0157/Claire_Stoll/RNASeq/11R-5_DLPRMzP300B4AAAPEI-2_2.fq

CONVERTING SAM TO BAM:
cd /fs/project/PMIU0157/Claire_Stoll/RNASeq/B73_hisat2/B73_hisat2_11-5
module load samtools
samtools view -bS a_B73_11-5_accepted_hits.sam > a_B73_11-5_accepted_hits.bam
samtools sort a_B73_11-5_accepted_hits.bam -o a_B73_11-5_accepted_hits.sorted.bam

cd /fs/project/PMIU0157/Claire_Stoll/RNASeq/Mo17_hisat2/Mo17_hisat2_11-5
module load samtools
samtools view -bS a_Mo17_11-5_accepted_hits.sam > a_Mo17_11-5_accepted_hits.bam
samtools sort a_Mo17_11-5_accepted_hits.bam -o a_Mo17_11-5_accepted_hits.sorted.bam

CUFFDIFF JOB
cd /fs/project/PMIU0157/Claire_Stoll/RNASeq/B73_hisat2/B73_hisat2_11-5
module load cufflinks/2.2.1
cuffdiff -b /fs/project/PMIU0157/Claire_Stoll/RNASeq/B73_hisat2/B73_hisat2_11-5 -u -N -p 8 /fs/project/PMIU0157/Claire_Stoll/Zea_mays.AGPv4.36.chr.gff3 /fs/project/PMIU0157/Claire_Stoll/RNASeq/B73_hisat2/B73_hisat2_11-5/a_B73_11-5_accepted_hits.sorted.bam /fs/project/PMIU0157/Claire_Stoll/RNASeq/Mo17_hisat2/Mo17_hisat2_11-5/a_Mo17_11-5_accepted_hits.sorted.bam

Is it the last command that gives the error?

I think cuffdiff is not quite giving you the right error...it looks like the crux of your problem is you aligned each sample to a different reference.

As @swbarnes2 noted, the error message is because the chromosome names, sizes, and orders differ between the two reference genomes/species you aligned against. Even if they were the same, you're comparing one sample against another between organisms, which is innately not something likely to work well. Further, I wouldn't expect the Zea Mays annotation to be very useful with whatever Ref_Lai is.


text:
Regular expression struggle
These are my files
SO_8543_IP_3_ablt_R1.fq.gz
SO_8543_IP_4_ablt_R1.fq.gz
SO_8543_INPUT_1_ablt_R1.fq.gz  SO_8543_IP_5_ablt_R1.fq.gz
SO_8543_INPUT_2_ablt_R1.fq.gz  SO_8543_IP_6_ablt_R1.fq.gz
SO_8543_INPUT_3_ablt_R1.fq.gz  SO_8543_IP_7_ablt_R1.fq.gz
SO_8543_INPUT_4_ablt_R1.fq.gz  SO_8543_IP_8_ablt_R1.fq.gz
SO_8543_IP_1_ablt_R1.fq.gz
SO_8543_IP_2_ablt_R1.fq.gz

Im running the below line
for i in $(ls *.fq*.gz | cut -f1-5 -d"_" | uniq);do echo $i; done

bowtie2 -q -p20 -x /run/media/punit/data1/BowtieIndex/hg38  -1 ${i}_R1.fq.gz | samtools view -bS - >${i%}.bam;

done

So if i use run this
for i in $(ls *.fq*.gz | cut -f1-5 -d"_" | uniq);do echo $i; done

I get files as such
SO_8543_INPUT_1_ablt
SO_8543_INPUT_2_ablt
SO_8543_INPUT_3_ablt
SO_8543_INPUT_4_ablt
SO_8543_IP_1_ablt
SO_8543_IP_2_ablt
SO_8543_IP_3_ablt
SO_8543_IP_4_ablt
SO_8543_IP_5_ablt
SO_8543_IP_6_ablt
SO_8543_IP_7_ablt
SO_8543_IP_8_ablt

The part of the file name that is consistent is _R1.fq.gz which goes after the input argument. As these are single end files so only 1 im giving one end as input.
Error message
***
Error: Must specify at least one read input with -U/-1/-2
(ERR): bowtie2-align exited with value 1
./Chip_SEQ.sh: line 25: syntax error near unexpected token `done'
./Chip_SEQ.sh: line 25: `done'

it looks like its not able to recognize as input .Any help or suggestion would be highly appreciated

You have an extra done, instead you want:
for i in $(ls *.fq*.gz | cut -f1-5 -d"_" | uniq); do
    echo $i
    bowtie2 -q -p20 -x /run/media/punit/data1/BowtieIndex/hg38 -1 ${i}_R1.fq.gz | samtools view -bS - >${i}.bam
done

As an aside, I'm not a big fan of piping ls to cut and then uniq. Why not instead:
for f in *.fq.gz; do
    sampleName=${f%%_R1.fq.gz}
    echo $sampleName
    bowtie2 -q ...
done

That will do the same thing, but (A) it's clearer what's going on and (B) it can be reused if your sample names change. Also, stop writing bash scripts and use snakemake or something like that instead, it scales much better and is less error-prone.

now after looking at my script number of times i figured out the "done" is extra .let me try your code. My shell scripting is very rudimentary so i use which i know is safe to save my self from screwing up .

i ran your first code the but i still get this error "Error: Must specify at least one read input with -U/-1/-2 (ERR): bowtie2-align exited with value 1" its  not able to recognize the input file, not sure what's wrong .

Try -U instead of -1 in the bowtie2 command.

yes that was the mistake i was doing thank you very much

@krushnachChandra please read this link to understand why you should basically never try to parse the output of ls.

@terdon I'd not say "never parse ls", I'd say "Don't parse ls unless you know exactly what you're doing, especially the filenames involved". I do a lot of for f in $(/bin/ls -1 <glob> | downstream_operations) because I'm 100% certain that each value of $f plays well with the other parts in the code.

@RamRS well, parsing ls is asking for trouble and just unnecessary. for f in <glob> does exactly the same thing but without calling an external process and without running any risk of failure because of file names. So why use a riskier and less efficient approach when it doesn't provide any benefit at all?

It's really useful when the $() feeding the loop involves a uniq, for example. It can't always be avoided. Whenever possible, use a glob. But don't have a hard constraint against $(/bin/ls | downstream_stuff)

thats a life saver link


text:
Normalization of data with RPkM
I'm having difficulty normalizing my data.
I was searching for transposable elements in my genome, and after this step, I made counts of reads in some transcripts. I produced something like this:
head(table_tissues_filtered_TE)
                  Lengths ova testes lobe retina suckers brain1 brain2 skin stage
Simple_repeat_80      134  58     77   48     69     115     137  131  195     75
tRNA_1                 59   0     14   12      1      19      12   14   21    104
Simple_repeat_87       26   1     33   12      3      15      24   21   19    180
Simple_repeat_114      22   0      0    0      1       0       0    0    2      7
Simple_repeat_115      30   0      0    0      0       0       0    0    0      1
Simple_repeat_123      22   2      3   317     45      13    652  651   15     21
                              axial                    gland viscera
Simple_repeat_80                 99                       35     557
tRNA_1                            9                        0       3
Simple_repeat_87                  9                        0       4
Simple_repeat_114                 0                      204       0
Simple_repeat_115                 0                       42       0
Simple_repeat_123               333                        5       4

where Lengths are the Length of each elements (simple repeats, etc), and the other columns indicate the reads counted with FeatureCounts.
I have another table with the number of reads for each tissue:
head(reads_table)
       ova   testes lobe   retina  suckers   brain1  brain2  skin   stage   axial
      522444 310243 226146  102307  126055   489389  668243  372728 262536  233754
  gland  viscera
  24817   25689

I would like to do an RPkM analysis to normalize the data using R, but I don't know exactly how to do it. Can anyone can help me? thank you!!!

It's very unlikely that "a RPKM analysis" is the right answer. Assuming you'd like to do differential expression, using tools like DESeq or EdgeR on the count table are likely to be a better thing to do.
For reasons why RPKM is not a good approach, have a read of this answer.


text:
Gene ratio as imput in limma
I have a data frame with gene-expression ratio. Is it possible to input this into limma/voom to find signinficannt gene-ratios between groups of samples?
my data:
                                   S1         S2         S3         S4         S5
hsa-let-7a-5p_hsa-let-7b-3p 80.066116 197.808511 275.322581 220.421053 242.100000
hsa-let-7a-5p_hsa-let-7b-5p  2.083889   2.393050   2.337716   1.062944   1.704930
hsa-let-7a-5p_hsa-let-7c-5p 83.517241  57.388889  85.350000  40.660194  71.205882
hsa-let-7a-5p_hsa-let-7d-3p  4.199393   5.792523   4.569058   3.193290   6.022388
hsa-let-7a-5p_hsa-let-7d-5p 32.729730  17.027473  14.490662   9.507378  22.416667

What do these ratio mean? Between what are those ratios calculated? But I'm not aware of any method that used ratios to calculate differentially expressed genes, as this is like a ratio


text:
SnpEff gene annotation with @ such as PCDHB@
I'm using snpEff to annotate some VCFs and getting some spurious annotations that don't have transcript ids. Some of these annotations have a gene name followed by ampersand such as PCDHB@ and I'm not finding any explanation in the docs.
Any idea what this refers to? Entire gene family?
Thanks.

The @ at the end of a name denotes a "gene cluster in chromosome region". This is mentioned in the official HUGO nomenclature guidelines in the "Characters reserved for specific usage" section:

It's also mentioned in the GeneCards documentation:

Very informative. Thanks!