How do I install hammock? I'm working with phage peptides libraries. Currently I'm analysing a massive data set from phage selection, NGS readout, preprocessed already (fasta file). I want to analyse the sequences to cluster peptides according to their consensus motifs. I found that Hammock may be helpful to do this.
Here is the GitHub repo for Hammock: https://github.com/krejciadam/hammock. I downloaded the master, unpacked it, but I don't know how to use it.
I guess I need to compile, but I have no idea what files need to be compiled, how to do it, what software I need use to compile etc.
I work with Mac/terminal.

Download a release (Hammock_v_1.2.0.7z) from https://github.com/krejciadam/hammock/releases
Unzip
To run the first example:

cd /your_path/Hammock_v_1.2.0/examples/MUSI/
java -jar ../../dist/Hammock.jar full -i musi.fa

As mentioned in the comments, it's in the manual and any further issues need more information (errors etc).


Run kallisto iteratively across many samples I am on a Mac using UNIX. I am trying to use the kallisto quant command on all files in a directory (instead of manually entering them). Because I'm running the analysis against the same index file, I first enter the following:

./kallisto index -i --index --make-unique index.fa

This successfully creates an index file. Then, I tried this:

for file in *.fasta; do kallisto quant --single -l 200 -s 0.1 -o $file-aligned
"$file"; done

While the $file-aligned folders are created, they're empty. I get an error that states that the index file is missing. So I assumed I would need to specify the index file in the command line. I then tried this:

for file in *.fasta; do kallisto quant -i index.fa --single -l 200 -s 0.1 -o
$file-aligned  "$file"; done

I get an error that says "Error: incompatible indices. Found version 0, expected version 10.

I have now also tried the following:

for file in *.fasta; do ./kallisto quant -i index --single - 200 -s 0.01 -o
"${file}-aligned" "${file}";

but I get an error saying the index file is not found.

I am now running

or file in *.fasta; do ./kallisto quant -i index --single -l 200 -s 0.1 -o $file-aligned

"$file"; done

This sort of works. The index file error messages are gone, but it runs quant on the same file 10 times instead of each file once. If I run echo $file it spits back the file that my command is running repeatedly.

If I close terminal, start a new session and run the same command, I get

Error: file not found ""

Running the jobs in parallel

Iteratively is one solution, but since you're not using bootstrapping, which can be run multithread (and is not needed if you're taking the popular tximport route) you can do better using GNU parallel, this will enable you to run as many jobs as you have execution threads on your system simultaneously.  Note that Kalliso only runs multithreaded with -t if you're using bootstrapping.

You'll need to make the index first:

kallisto index --make-unique -I index.fa.idx index.fa

Next lets define the input FASTQ files and the sample IDs you want to use using a tab-delimited input file like so:

sample_1    sample_1_R1.fastq.gz
sample_2    sample_2_R1.fastq.gz
sample_3    sample_3_R1.fastq.gz

Here the first column is the sample ID which will be used for output and the second separated by a tab character is the input file.

You can produce this using a text editor but as a short-cut you can feed the list of fastq.gz you have into a list via:

ls -1 *.fastq.gz > editme

Next we create a list of jobs to run on our system using a utility script Make_job_list_kallisto.sh which consumes our tab delimited input file list above.  I've used your -l and -s parameters below, note that with paired end data normally these are estimated for you.

#!/bin/bash -eu

[ $# -ne 2 ] && { echo -en "\n*** This script generates jobs for GNU parallel. *** \n\n Error Nothing to do, usage: < input tab delimited list > < output run list file >\n\n" ; exit 1; }
set -o pipefail

# Get command-line args
INPUT_LIST=$1
OUTPUT=$2

# Set counter
COUNT=1
END=$(wc -l $INPUT_LIST | awk '{print $1}')

echo " "
echo " * Input file is: $INPUT_LIST"
echo " * Number of runs: $END"
echo " * Output job list for GNU parallel saved to: $OUTPUT"
echo " "

# Main bit of command-line for job
CMD="kallisto quant -i index.fa.idx --single -l 200 -s 0.1"

# Main Loop
[ -e $OUTPUT ] && rm $OUTPUT
while [ $COUNT -le $END ];
do
    LINE=( $(awk "NR==$COUNT" $INPUT_LIST) )
    # Make file list
    echo "Working on $COUNT of $END Sample ID: ${LINE[0]}, Files ${LINE[@]:1}"
    echo "$CMD -o ${LINE[0]} ${LINE[@]:1}" >> $OUTPUT
    ((COUNT++))
done

We use the script like so:

 Make_job_list_kallisto.sh input_file_list.txt job_list.txt

Now we have a job list for GNU parallel we can run our jobs in parallel like so:

parallel --progress --jobs 4 --joblog kallisto_joblog.txt < job_list.txt

This will run 4 of the kallisto jobs simultaneously (Increase this if you have more threads - don't mind consuming all the IO/CPU threads on your system) producing a nice progress report of those completed, additionally a job log will be written to the file kallisto_joblog.txt.

For paired end data I already have these scripts + similar instructions on GitHub which people might find useful.


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

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

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

plot

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

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


Microbial Classification with 16S Nanopore Data I have recently been tasked with analyzing a heterogeneous microbial sample composed of 16S data generated on a MinION sequencer, with the goal of assigning taxonomies and calculating respective abundance.

It is my understanding that traditional preprocessing methods such as Dada2, Unoise3, and Deblur would likely not be appropriate for this kind of data given the poorer read quality inherent with nanopore sequencing. To give an example, the average read quality of this current experiment is a Phred score of 15 and an average read length of 1500bp.

Furthermore, it would seem best to try and utilize a classifier that makes use of a full length 16S reference, rather than just the V3 region, to try and offset the higher error rate. But are there any classifiers that do this out of the box?

I am very interested to hear of any advice or input the community has on this. Thanks!

Use kraken2 with one of the 16S databases. Kraken2 has worked well for me with nanopore sequences. Pavian can be used to inspect the Kraken2 report files visually.

Using uncorrected called reads has worked for me. Any errors from nanopore sequencing should be inconsistent enough (i.e. not always lead towards the same wrong organism) that the noise doesn't matter.

Note: if you're using the MiniKraken database instead of a comprehensive 16S database, there might be some false positive species assignments depending on which kmers get pruned out.

Another possibility is centrifuge, but I'm not sure if that works for 16S amplicon sequencing; it's primarily for metagenomic analysis.


What algorithms exist for reverse protein folding? Are there any algorithms already in existence for reverse protein folding/docking (e.g. taking a desired surface map/docking profile and generating DNA which will produce a protein to fold into it?)

Protein folding seems straightforward (though computationally intensive) and protein docking is likewise straightforward, but thus far I haven't come across any protein unfolding algorithms, though I might just be searching for the wrong keywords.

Any advice in this area?

The researchers behind FoldIt have the occasional design challenges that use human game players to help design proteins to fit a structure.

However, even with a known protein sequence there's still a lot of labwork involved in trying to find a good working DNA sequence. Just generating the DNA sequence is a difficult task due to redundancy in the genetic code (i.e. 64 codons for 20 amino acids + stop signals). When that is combined with cell-specific tRNA frequencies that can influence the timing of amino acid incorporation (and therefore how a protein folds as it is created), creating a sequence that reliably fits a shape ends up being quite difficult.


classifying samples by TCGA signature I have some RNA-seq samples from multiple glioblastoma tumours that I'm now trying to classify according to a specific gene signature (from Verhaak et al., 2010) using R. The gene signature is reported as a gene list with specific centroids for each of the 4 clusters (https://api.gdc.cancer.gov/data/941f81a1-05d7-4f84-80ec-534b8dc1ebac). I'm wondering how I can use this signature to classify my samples in R? Would it involve some sort of k-nearest neighbours method?
Additionally, the signature was identified using microarray data, but I am classifying RNA-seq data. Is there any sort of adjustment I should make to the signature to account for this?
Thanks in advance!

This is not a very simple problem. I saved the table and "test.txt" and visualize it below:
library(pheatmap)
x = read.delim("test.txt")

head(x)
  Gene.Symbol            Composite_chr_coords  Proneural      Neural
1      CDKN1B   [36.1:12:12761576-12766570:+] 0.03353670  0.03353670
2       EPB41    [36.1:1:29086215-29319545:+] 0.03882308  0.03882308
3        CLGN  [36.1:4:141529057-141568232:-] 2.12796056  0.41999847
4      PDE10A  [36.1:6:165664525-165995575:-] 1.66012916  0.25349464
5     RALGPS2  [36.1:1:176960923-177152465:+] 0.78971500 -0.01945154
6        TAF5 [36.1:10:105117714-105138812:+] 1.08354870  0.09461930
    Classical Mesenchymal
1  0.03353670 -0.84379400
2 -0.84821990  0.03882308
3  0.41999847  0.41999847
4  0.25349464  0.25349464
5 -0.01945154 -0.01945154
6  0.09461930  0.09461930

You can see there are no so called markers genes and if you look at the scale of the data, it was normalized in microarray and it will be very hard for you to get your RNA-seq to fall into that range, and apply the medoids to identify your clusters or type.
Since you did not provide any data, and it is not trivial to simulate data for this, I can suggest the approach in words here:

Subset your count table to these exact set of 840 genes.

Normalize the data and do a PCA or hierarchical clustering to see if you indeed observe 4 or 3 or 2 subtypes

If you indeed see 4 or 3 clusters this would be straight forward, take the centres of means of these clusters, and correlate with each of the sub types here. Chances are the subtype that is most correlated with one cluster, is the subtype you are interested in.


How do I force Cromwell / WDL to put input files in the same directory? I have a cromwell task that, among other variables, takes these three variables:
   File Ref             # Reference genome fasta
   File RefFai          # Reference Genome index
   File RefDict         # Reference Dictionary

While all three files make it into the task's "inputs" directory, they end up in different sub-directories of inputs:
ME $ find inputs
inputs

inputs/-181926642
inputs/-181926642/Sample.vcf.bgz.tbi
inputs/-181926642/Sample.vcf.bgz
inputs/-1548908437
inputs/-1548908437/GCA_000001405.15_GRCh38_no_alt_short_headers_nonACTG_to_N.dict
inputs/44588251
inputs/44588251/GCA_000001405.15_GRCh38_no_alt_short_headers_nonACTG_to_N.fa
inputs/44588251/GCA_000001405.15_GRCh38_no_alt_short_headers_nonACTG_to_N.fa.fai

Which means that when I call GATK and pass it the path to the .fa file, it gets rather cranky that the .dict file isn't in the same directory.
How do I force cromwell to put all three files in the same directory?

My hacktastic fix for this is to use basename to get the get the name of each file, then do
refName=`basename "${Ref}"`
ln -s "${Ref}" "$refName"

for each file.  This puts aliases to each file in the executions dir, and then I pass "refName" to GATK rather than "${Ref}".
Not a good solution, but a solution


How to QC overamplified shotgun library? I made some NEB ultra II whole-genome shotgun libraries with around ~200ng of template going into the indexing PCR step. The PCR was run for 13 cycles and I ended up having high library concentrations between 60-90 ng/µL. Should I be worried about complexity or PCR duplicates? How would all of you QC to make sure the library is complex enough and not overamplified too badly?

FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) is a good start: it will show you observed/expected ratios for the most overrepresented sequences, and running it on a few datasets you know to be high quality should give you intuition for what looks good and what doesn't.

My labmate also produced something called preseq (http://smithlabresearch.org/software/preseq/) that uses the library complexity of an initial sequencing experiment to predict how many more unique reads you would get with additional sequencing of the same library. It could give you some intuition for whether or not it is worth sequencing your library deeper to get the information you need.


Do I need to remove water molecules before protein-ligand docking? I am going to screen small molecules against a protein target. Do I need to remove all the solvent molecules before performing VLS?

Yes, if you are using implicit solvent you ought to remove all of them.
Exception. Unless the crystallographic waters are trapped essential structural components of your active site and your ligand is not going to replace it, in which case keep it because its remove changes the energy of the active site —say you have a catalytic triad, with a water that is pinned between the acid (e.g. aspartate) and the base histidine, its removal would change the tautomerisation of latter and result in the nucleophile being protonated (and not nucleophilic). Having said that, if you don't have time to look at your structure (or cannot because of automation) and are okay with some artefacts from incorrect protonation strip them.
If you are not using implicit solvent, then it depends on the tool used —most require them stripped and re-added. If you do keep them, make sure that the crystallographic water (HOH) has the same residue code as the added water (TP3, WAT etc.).


ncbi...eutils/esearch.fcgi: How to search now in the snp database? According to the 2017(!) textbook, it worked at that time, for instance,

https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=snp&term=autism

At this time, this does not work and for each conceivable term it will find 0 records.

What is the right way to search in snp, valid in December 2018?

eUtilities for dbSNP is in the process of being retired. See this post from NCBI. You may want to check out the new API.
From the NCBI Insights blog post:

Entrez search is available only for human. We will retire all Eutils report formats (RS docsum XML and ASN.1, FASTA, and FLAT text) later this year. Only eSummary XML report will be supported.
Batch Query Service for all organisms will be retired October 2018. Users of human dbSNP data should transition to the new API.
We will place non-human FTP data in the archive directory. Please see the phase-out plan announced last year.


What's a good ontology for drug names? I have ...

A database with patient phenotypes in, stored as HPO terms
Genetic data in whatever format I need

I want ...

To store drug names in a way that won't make my life difficult
If it's feasible, to make them so that you could draw links between HPO terms, medications, genetic changes, and outcomes.  Note that I'm not actively trying to do that, but I don't want to exclude the possibility of it in the future

Ideally, I'd like those drug names in some kind of ontology, and even more ideally an ontology that's available on http://bioportal.bioontology.org

So - what's the best drug ontology to use in this case?

Update

  What relationships would this drug ontology describe? Structural similarity of the compounds? Affected pathways? How would they be linked? Do you want all drugs that can cause the same side effect linked, for instance? Or all drugs that can be used for the same disease?
  Drugs by disease would be the most useful.  The information that would be most useful would be something like

Atenolol

Atenolol 25mg
Atenolol 50mg
Atenolol 100mg

Propanolol

Propanolol 10mg
Propanolol 40mg
Propanolol 80mg
...

In short, being able to get a list of dosages for a given drug.  It'd be extra-useful if I could also retrieve a set of drugs for a given condition, so in the above example to retrieve a list of beta blockers (Atenolol, Metoprolol, Nebivolol, ...)

  Have you looked at any of the various drug ontologies listed in bioportal.bioontology.org? Which ones have you discarded so we don't suggest the same ones? What were they missing that you need?

MDDB doesn't seem to be structured as above
Likewise The Drug Ontology
RxNorm seems to be the best available that I've seen

Since this question's been asked there's been a (what looks like a first useful?) release of Drug Ontology DRON. Perhaps worth checking again?


What does presence/absence of black lines mean in IGV? IGV (Integrated Genome Viewer) is a popular open-source tool for viewing alignment files. In my BAM file in IGV, some deletions have black lines running through them, and others don't. What causes this? Short black lines also appear adjacent to some mutations and insertions.
I generated the CIGAR for the SAM file myself, and am wondering if I made an error. I didn't find anything in IGV's documentation about this.

Black line represents an indel.  Read was lengthened to fit consensus.
No black line (blank), no read/sequence.  Read likely ended or was truncated.


how to remove range from fasta header Could you please suggest me how I can remove range from fasta header:like these number from below sequences which has some colons indicating the range of the genome

:147010-147657

:149201-149845

> FVJT01000004.1:147010-147657 Mycobacterium abscessus subsp. massiliense strain 974
GTGACCAGCAGCGCAGACCGATCCCCCGTCGAGGCGATGGTGAGCCACGCCGAGGGCGCCATCTCCGAGG
ATGACATCGTCGCCGCCGCCCGGGAACGTGCCGTTGACCTGGGCGCGGGCCCGGTGACACCCGCGGTAGG

>FSAX01000002.1:149201-149845 Mycobacterium abscessus subsp. bolletii strain 286 genome
ATGACGACGCCCGTCAGCTTCGATAACAGATTCTTCTCCGCAACGTGGTCTGCTATCGCGCGGCGCGAAC
CTCAGGAGATCAAGGAGCTGCGTGCGGAGAATCTACGCGGCCTGACCGGCCACGTTCTGGAGGTGGGAGC

Invest some time in learning unix tools such as grep, sed and awk.

sed -r '/^>/s/:[0-9]+-[0-9]+//' input.fasta > output.fasta

This line of code matches the pattern :<numbers>-<numbers> on lines that start with > (FASTA sequence headers) and replaces that pattern with '', essentially removing text that matches the pattern.


How to build a BAM header file with htslib in C++? I'd like to use C++ to generate a new BAM file programmatically. This is an example how to use htslib to generate a new BCF file on the fly.

  https://github.com/samtools/htslib/blob/develop/test/test-bcf-translate.c

The two key functions are:

bcf_hdr_append
bcf_hdr_add_sample

Now when I navigate to the header file for SAM/BAM functions:

  https://github.com/samtools/htslib/blob/develop/htslib/sam.h

The only "append" function is bam_aux_append. But that's for optional fields.

My questions:

Given a bam_hdr_t header (potentially empty), how to add a new chromosome reference? (For example, add a new chrAB to the BAM file)
Given a non-empty bam_hdr_t * header, how to remove a chromosome? (For example, remove chr4 from the BAM header)

My sample code:

// Empty BAM header, nothing inside
bam_hdr_t * hdr = bam_hdr_init();

// I'd like to add a new synthetic `chrAB` into the header, how?

You're just looking in the wrong header file. The BCF stuff will be in vcf.h (you'll see the two functions you mentioned there).
With BAM files one doesn't typically need to add/remove headers, which is why there are not convenience functions for them. To add a header:

Increment hdr->ntargets
realloc(hdr->target_len, ...) and hdr->target_name.
You likely want to append the name/length to the end of hdr->target_name and hdr->target_len, respectively.

I don't think you need to update the dictionary or plain text for most uses (or at least I've never needed to).
Removing a chromosome would follow the reverse instructions, (you'd need to shift everything in the name and length arrays if you remove entries in the middle.
WARNING, doing this is prone to error!
I would strongly cautious you to be very careful when doing this. Please remember that the chromosome that entries are mapped to isn't actually stored in alignment entries themselves. Rather, they simply have an index into hdr->target_name, which saves space. Consequently, if you start adding/removing headers to/from a BAM file with alignments then you're going to have to modify all of the alignments to keep the correct alignment:chromosome associations. Otherwise you'll end up with completely useless results (this is the same caveat as with samtools reheader).


Can gene co-expression networks be used to help identify differentially expressed genes? In my RNAseq dataset of differentiated stem-cell lines, some samples have far fewer significantly differentially expressed genes than others. QC shows that this is because there are way fewer reads for one sample than the others.

Can gene co-expression networks be used to infer differential expression? For example, if certain genes just miss the cutoff q-value cutoff  If so, what papers have done so?

While you can use networks to find differentially expressed genes (see the WGCNA package, which does this) in my experience this ends up largely matching what you'd get using a traditional package with a looser threshold for significance. Given the time savings of traditional packages, there's rarely any gain to using networks (it won't hurt to try, just note that it'll take some time). If some very interesting genes are just above your significance threshold then change your threshold. You need to validate your findings in some way anyway, so your p-value threshold is partly just a way to protect you from wasting your time (but if background knowledge suggests that your time wouldn't be wasted...).


How to handle control samples in CLIP-Seq I have a CLIP-Seq dataset I'm processing, which includes control samples and no inputs.  This is the second CLIP analysis I've performed to help out users of our genome core facility and the first one only had inputs. I've done a number of ChIP analyses in the past, which I loosely consider to be a similar type of analysis, and they all used inputs, so I've never had to look into the differences in using inputs versus controls when doing peak calling before.  Now that I'm trying to figure out what I need to do, I'm a little lost.

Here's what I've got:

Samples of cross-linked RNA/protein complex that went through a poly-A and then a FLAG pull-down (protein G beads with an antibody).  The complex was eluted, treated with proteinase K, and the RNA was sequenced.
Controls went through the same processing as the samples except the control beads had no antibody

There are replicates in both cases.  (As a side-note, the RNA was treated with RNase III while bound to the beads and the RNA's trimmed off the complexes were also sequenced.  I'm not sure what to do with those samples yet.)

I was going to use one of the CLIP peak callers available on Galaxy, but whichever one I use, all the protocols I've found and the documentation I've looked at use inputs for the peak calling method and not controls.  So what do I need to do differently (if anything) to use the controls as the background in the peak calling?  Or do I need to suggest they also sequence the inputs?

And just so I know (dumb question here)... Is "input" (as the name suggests) the sample as it was before the FLAG pulldown?  Why would one use one or the other (input versus control).

I've not handled CLIP-seq before, but we have done iCLIP and in this we had controls, but no input. @DevonRyan is correct, these control samples are similar to IgG pulldowns in ChIP seq.

The way we handled this was to call our peaks on isolated samples. The original peak calling alog from the first CLIP seq paper doesn't require input. Also I believe that clipper doesn't need input (although it will take it). We then removed any peaks from the sample that overlapped a peak in the control.

But most of our analyses we didn't do on peaks - we did lots of things like looking at metagenes, or classes of genes, or particular types of exon etc, and when we did this, because we were averaging over many sites, we could generate the pattern for the control as well as for the input, and often normalise the one to the other.


Calculating read average length in a Fastq file with bioawk/awk I found here this awk script:

BEGIN {
  headertype="";
  }
{
  if($0 ~ "^@") {
    countread++;
    headertype="@";
  }
  else if($0 ~ "^+") {
    headertype="+";
  }
  else if(headertype="@") { # This is a nuc sequence
    len=length($0);
    if (len>4) {
      readlength[len]++;
    }
  }
}
END {
  for (i in readlength){
    countstored+=readlength[i];
    lensum+=readlength[i]*i;
    print i, readlength[i];
  }
  print "reads read = "countread > "/dev/stderr";
  print "reads stored = "countstored > "/dev/stderr";
  print "read average length = "lensum/countstored > "/dev/stderr";
}

and I just wonder if it is possible to shorten it with bioawk?

This can also be done with regular awk.

awk '{if(NR%4==2) {count++; bases += length} } END{print bases/count}' <fastq_file>

The NR%4==2 count the second line out of every block of 4.  length is a built-in that defaults to the length of the line, same as length($0).  In this case you can inject you custom printing to the END{} block but countread and countstored will always be the same since the averaging is done on the fly.


How to merge sam files together with adding read groups I have three sequencing libraries of single individual mapped to a reference using bwa-mem. I would like to merge the three unsorted .sam files I have so, I can call variants and heterozygosity estimates using atlas. Atlas requires one input mapping file (bam) with defined read groups because it estimates the error profiles of different libraries separately.

How can I merge multiple sam files? Preferably avoiding java (Picard tools).

I tried to figure out a solution using samtools 1.3. I sorted individual files using samtools sort, then I used samtools merge -r merged.bam s1.sort.sam s2.sort.sam s3.sort.sam to merge the sorted files. However, the read group didn't make it to the header (and the variant caller I use is complaining about it), also the read group is stupidly the file name.

I tried to define meaningful readgroup names using procedure described on BioStars, but I found that this will just change the header, it does not adjust the names of read groups defined by samtools merge (the file names).

Following this related thread on SeqAnswers, I tried to define the correct header with read groups corresponding to merged file names:

samtools -rh rg.txt merged.bam s1.sort.sam s2.sort.sam s3.sort.sam

where rg.txt is a file with header

@RG s1.sort
@RG s2.sort
@RG s3.sort
...
output of samtools view -H s1.sort

However, the header still had not the read group, I guess because sam header accepts only tagged items to be specified (something like @RG    XY:s1.sort). So, I looked into the merged bam and I found that the tag of RG is Z:. So I tried to just rename header of the merged file using samtools reheader, but then samtools complain about the fact that a tag needs to be of length 2:

Malformed key:value pair at line 123: "@RG    Z:s1.sort"
Segmentation fault (core dumped)

I have opened an issue to report this strange incompatibility of read groups generated by samtools merge with samtools reheader.

I would like solution to :

create a bit more standardized read group names (SM:Sample\tLB:library format)
avoid pointless writing to disk like in sam -> sorted.sam -> merged.bam case (can be probably achieved through "pipes and tees", thanks @bli)

I also know that I can specify RG to bwa, so the sam files will have read groups defined in the first place. But I do not like the idea of remapping three libraries just to create correct formatting of read groups.

I have used https://github.com/ekg/bamaddrg to add read groups quickly to multiple sam files.  And then you can do a samtools merge of the tagged files.


Can open-source software be peer-reviewed and published? My colleague and I have developed a software tool intend to release it open-source. This tool is specifically for tasks in bioinformatics but we think it would be helpful for the wider community. Our institution will permit us to release it provided we get appropriate credit.

Thus we wish to publish it in peer-review. Is peer-review publication of bioinformatics software available? If so what is required to publish it?

There are whole journals based primarily around publishing open source tools. The primary example of that is "Bioinformatics", where a lot of the open-source tools are published. We've also had luck publishing in the Nucleic Acids Research yearly special webservers issue, since we make Galaxy wrappers around our tools.

You can also publish on tools in venues like Nature Methods, but it's vastly more difficult to get in there if your paper is purely on a tool and it's not something game-changing like salmon or kallisto.

Regarding what's actually needed to publish a tool, you will typically need the following:

Decent documentation. I ask for this as a reviewer and am asked about it when I submit things.
Easy installation. No one wants to spend two hours trying to compile your tool, so make sure there's a conda package or at least a docker container for it.
Test data. Without this your reviewers and users will likely not be happy since they won't be able to try it out on something small.

For the actual submission you'll tend to need comparisons to other tools. I'm generally not a big fan of this, but you will find that some reviewers will demand it (assuming there are other tools that do something similarish to yours).


After running a DRIMSeq pipeline, how do I know which genes are upregulated in the different conditions? After running a DRIMSeq pipeline and obtaining the genes that are differently used between the null and full models, how do I select the genes that are differently used in the different conditions?

According to the vignette, I would have to use the plotProportions function and say in the group_variable the variable of my interest. What if I want a list of the genes that are differently used specifically in the variables of my interest? Checking every gene manually sounds exhausting. Help would be much appreciated.

After running the test for differential transcript usage you get a list with gene- and transcript-level p-values that you can use to filter which differences are significant. The gene-level p-value tells you if there is differential usage of transcripts within a gene and the transcript-level p-value tells you wich transcripts actually are changing.  Therefor you will want to select significant genes for further investigation.

The plots are mainly used to visualize top hits as detailed in the vignette.

res <- results(d)
# The data is ordered by p-value
res <- res[order(res$pvalue, decreasing = FALSE), ]
# The top hit is selected
top_gene_id <- res$gene_id[1]
plotProportions(d, gene_id = top_gene_id, group_variable = "group")

However you should consider running the two-stage test that is described on page 16 of the vignette you are using to reduce the number of false positive hits. From this you will get another list of adjusted p-values that can be used as above.


DMPfold: number of models produced using "10 500" feature in "run_dmpfold.sh"? I ran this command:

run_dmpfold.sh sequence.fasta sequence.21c sequence.map ./folder 10 500

From my previous discussion, I understand that I should run with default "run_dmpfold.sh" parameters i.e. 3 50; however, since I have already run the 10 500 I have some questions.

I see that after generating "sequence_1.pdb to sequence_500.pdb" files, these files are no more in the folder and, now I see "ensemble_1.pdb".
Will there be more ensemble files like "ensemble_500.pdb" and "sequence_n.pdb" files? I mean when can I expect "final_1.pdb' file? Also, I see files like "CLUSTERS_001.pdb".

The sequence_n.pdb files are temporary and are combined into an ensemble.i.pdb file after each iteration i. That file has many structures in. After the run is complete there will be a final_1.pdb file. Don't worry about the cluster files. If you are interested in what they do, consult the run_dmpfold.sh script and see what it is doing line by line.

As I have said multiple times elsewhere, you should run the script with default parameters and let it run to completion. Then use the final_1.pdb and any other final_*.pdb files as the output. If you want more structures, use the ensemble.i.pdb files.


>My counter is counting genotypic combination occurences more than once, how do I ensure it counts one combination and doesnt count it again? mutations = ['A222V', 'D614G', 'E484Q', 'E780Q', 'G476S', 'L18F', 'N439K',
             'S477', 'S477N', 'T478I', 'V483A']

combinations = []
for M in range(1, len(mutations)+1):
    for subset in itertools.combinations(mutations, M):
        combinations.append(subset)

        combinations = ['_'.join(sorted(x)) for x in combinations]
        combinations = [x.split('_') for x in list(set(combinations))]

root = "C:"

os.chdir(root)

lineages = os.listdir('Results')

combination_labels = []
combination_counts = []

for lineage in lineages:

    df = pd.read_csv('Results/' + lineage).dropna()

    for combination in combinations:

            combination_df = df[list(df)]

            for mutation in combination:
                combination_df = combination_df[combination_df[mutation] == 1]

        #print(combination_df.shape[0])

            combination_labels.append('_'.join(combination))
            combination_counts.append(combination_df.shape[0])

            out_df = pd.DataFrame({'combination':combination_labels,
                               'count':combination_counts})
            out_df['percentage'] = (out_df['count'] / df.shape[0]) * 100
            out_df = out_df.sort_values('percentage', ascending = False)

            out_df.to_csv('Results_2/' + lineage.replace(".csv", "") + '_3.csv',
                     header = True,
                     index = False)

The input CSV
lineage,A222V,D614G,E484Q,E780Q,G476S,L18F,N439K,S477,S477N,T478I,V483A
417941,0,1,0,0,0,0,0,0,0,0,0

Output CSV
combination,count,percentage
D614G,87355,90.7084929856806

My above code is used to count all occurences of combinations of spike protein mutations

to determine the genotypes.
However, looking at the values processed into the .csv file, i can see ' D614G' is being counted even when
it is in combination with other combinations.

My question is how do i ensure once something has been counted (a row in the imported .csv file)
it is not counted for a second time?
OR possibly how can I edit this code to prevent counting of singular mutations even when presented within a combination?
Thanks in advance

Alright, so there are a number of problematic patterns in your code - as far as I understand what you are trying to do. Next time, try to post a reproducible example that people can use and more people will be willing to help.
L18-21:
combination_labels = []
combination_counts = []

for lineage in lineages:

Declaring these two lists before the loop, then appending all your combination labels inside the loop means that these two lists will contain all labels and counts for all runs (all your lineages). If I am understanding your intentions correctly this is probably not something you want. The simple fix here is to just move them into the for-loop.
L29-30:
for mutation in combination:
    combination_df = combination_df[combination_df[mutation] == 1]

Here you are iteratively checking if the value for each mutation is set to '1', but you are also overwriting your combination_df variable at each iteration of the loop. After this loop, combination_df will be whatever mutation last was set to '1' in your current combination. I will get back to solving this in my solution below.
L37-44:
combination_labels.append('_'.join(combination))
combination_counts.append(combination_df.shape[0])

out_df = pd.DataFrame({'combination':combination_labels,
                   'count':combination_counts})
out_df['percentage'] = (out_df['count'] / df.shape[0]) * 100
out_df = out_df.sort_values('percentage', ascending = False)

out_df.to_csv('Results_2/' + lineage.replace(".csv", "") + '_3.csv',
         header = True,
         index = False)

Since this part is again scoped inside the for-loop, each iteration (each combination) will overwrite the same 'Results_2/' + lineage.replace(".csv", "") + '_3.csv' file. This you also want to move out of the for combination in ... loop.

A big concern is also performance. You have 2000-something combinations and 100,000 lines in your file. As you loop over combinations, you also process the DataFrame each time, when this is a proble you can solve in one iteration. It's not much data so it's still fine, but these habits are still good to develop early.
Now here's an idea for a rewrite that should work better. I don't have a csv from you to test with, so there might (proabably are) still things to fix, but hopefully it'll give you a starting point:
mutations = ['A222V', 'D614G', 'E484Q', 'E780Q', 'G476S', 'L18F', 'N439K',
             'S477', 'S477N', 'T478I', 'V483A']

combinations = []
for M in range(1, len(mutations)+1):
    for subset in itertools.combinations(mutations, M):
        combinations.append(subset)

combinations = ['_'.join(sorted(x)) for x in combinations]

root = "C:"

os.chdir(root)

lineages = os.listdir('Results')

for lineage in lineages:

    out_df = pd.DataFrame({"Count": 0}, index=combinations)

    df = pd.read_csv('Results/' + lineage, index_col="lineage").dropna()

    for index, row in df.iterrows():
        combination = '_'.join(sorted(df.columns[row==1]))
        out_df.loc[combination, "Count"] += 1

    out_df['Percentage'] = (out_df['Count'] / df.shape[0]) * 100
    out_df = out_df.sort_values('Percentage', ascending = False)

    out_df.to_csv('Results_2/' + lineage.replace(".csv", "") + '_3.csv', header = True,
        index = True)


Clustered column chart for single cell data I want to create a column plot for my single cell analysis that has 2 sample types, a normal and PDAC. I want to show, the percentage of different cell types per sample type.

How do I go about determining the percent of each cell type in each condition and then go about plotting the column plot? Also, is it possible to determine the SD from the dataset and add it?

dput(head(all@meta.data))

structure(list(orig.ident = c("PDAC", "PDAC", "PDAC", "PDAC",
"PDAC", "PDAC"), nCount_RNA = c(7945, 7616, 7849, 5499, 853,
1039), nFeature_RNA = c(2497L, 2272L, 2303L, 2229L, 509L, 588L
), percent.mt = c(2.63058527375708, 3.7421218487395, 4.9433048796025,
0.490998363338789, 0.234466588511137, 0.192492781520693), RNA_snn_res.0.5 = structure(c(5L,
5L, 5L, 5L, 1L, 6L), .Label = c("0", "1", "2", "3", "4", "5",
"6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16",
"17", "18", "19"), class = "factor"), seurat_clusters = structure(c(5L,
5L, 5L, 5L, 1L, 6L), .Label = c("0", "1", "2", "3", "4", "5",
"6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16",
"17", "18", "19"), class = "factor"), S.Score = c(0.0171539298241524,
-0.00599651468836135, 0.011058574403656, -0.0424413740104151,
0.03549521163095, -0.0243347616753686), G2M.Score = c(-0.0591387992008088,
0.0253275642795205, -0.0512402869839816, -0.0239076248512967,
-0.0263164126515597, -0.0339086517371782), Phase = structure(c(3L,
2L, 3L, 1L, 3L, 1L), .Label = c("G1", "G2M", "S"), class = "factor"),
    old.ident = structure(c(1L, 1L, 1L, 1L, 1L, 3L), .Label = c("Fibroblast",
    "T Cell", "Endothelial", "Tumor", "Stellate", "Macrophage",
    "B Cell", "Mast Cell", "Acinar", "Endocrine", "Exocrine"), class = "factor")), row.names = c("PDAC_AAACCCAGTCGGCTAC-1",
"PDAC_AAACGAAGTCCAGGTC-1", "PDAC_AAAGAACAGCAAGTGC-1", "PDAC_AAAGAACAGGATGAGA-1",
"PDAC_AAAGTGACATCCGTTC-1", "PDAC_AACAACCAGGAAGTAG-1"), class = "data.frame")

You can make use of table() and prop.table() to calculate the number and proportions of your cell types with respect to tissue of origin.

prop.table() requires a table as input so you can start with something like prop.table(table(object@meta.data$cell_type, object@meta.data$tissue). This will create an objet of class table.

Next step is to create a data frame with as.data.frame(your_table_object) and then to convert it to "long format", for example with tidyr::pivot_longer(). You can feed these data in long format to ggplot().


GSEA preranked different results opposing phenotypes I am running an analysis of two phenotypes: mutant and control.
I have a statistic used to rank the different genes, and the ranking is exactly the opposite between mutant and control. Moreover, there are no duplicates in the ranking (which excludes that genes are ranked alphabetically). I uploaded the test file here.
My expectation is that running GSEA using mutant as positive phenotype, should give me the same results (but in opposite direction) than running GSEA using control as a positive phenotype - since there are only two phenotypes, and the ranking statistic is computed mutant vs control.
However, when I run Gene Set Enrichment Analysis using mutant as positive phenotype, the results I get are not the same with respect to running GSEA using mutant as negative phenotype.
I initially thought that it might be related to the choice of the random seed, however, I get the same results even setting the same random seed in all runs.
The following is an example of an enriched result of the positive phenotype:

The same results with the negative phenotype:

The ranking is exactly the same, however, the leading edge is different in the two cases. It's probably due to the fact that the ranking sum goes in one direction (increasing the enrichment score) in both phenotypes. This is probably the cause of different scores in the opposite phenotypes, however, which is the result I should use in this case?

What do you mean by having the same results but opposite direction?
What GSEA does is that it goes through the gene list from the top to bottom, whenever it encounters a gene that belongs to a gene set, it adds a positive score to that gene set, and if not, it penalizes the score. That is what you see in the green line in the results.
Now imagine that you reverse the order of your input gene list, it is obvious the results are going to be different if the program still goes through the list from top to bottom.
I think as a general rule, whenever you change the order of input gene list, results are going to be different.


Mining documents for nucleotide and peptide sequences The text mining literature has an emphasis on identifying and normalizing gene names, mutations, pathways, concepts, and so on. I haven't been able to find much, however, on methods for extracting nucleotide and peptide sequences from documents. I did find this tidbit from Aerts et al. (2008), emphasis mine:

Text was split into words and words greater than 10 characters in length with greater than 40% of characters from the capitalized DNA alphabet [ACGT] were extracted using regular expressions to isolate putative DNA sequences.
All putative DNA sequences extracted from each paper were concatenated in the order they appeared in the text into a single fasta sequence and labeled with the corresponding PMID.
Concatenation of sequences was performed to merge sequences split by line breaks in the text conversion, and because we reasoned that inappropriate joins would be reconciled at the genome level by local alignment procedures.
Extracted, concatenated sequences were used as queries to BLAST RepeatMasked versions of genome sequences downloaded from the UCSC genome database...

I'm struck by how simple the author's method for extracting DNA sequences is. But I guess with the limited complexity of the DNA alphabet, that makes sense. This approach wouldn't work with the much more complex protein alphabet.
Are simple regular expressions the text mining community's state of the art for finding DNA/RNA sequences in documents? What methods are used to identify protein sequences?

I am not sure there is a state of the art for that.
Some possible strategies I can think of:

Applying a dictionary to identify language words and discard these.
Using some context information like identifying some words in the same sentence than the sequence.
I guess the format in papers for sequences will not be the same as for text. Maybe you can find a pattern.

I would use an incremental strategy. Use the simplest algorithm then see what it recovers. If there is garbage, find a strategy to remove some garbage with a common pattern. Repeat until you are happy with the result.
Disclaimer: I have some experience in Natural Language Processing.


Comparing RNA-Seq and Ribo-Seq I'm attempting to compare total RNA-seq with Ribo-Seq, to determine if changes in Ribo-Seq are due to changes in transcriptional expression (akin to the analysis performed by Anota2Seq). I am, however, struggling to understand how to accomplish this in DESeq2.

My design is as follows:

cts <- read.delim("Interaction.txt", row.names=1)
coldata <- read.delim("Sample_Sheet.txt", row.names=1)
cts <- round(cts)
dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~Genotype + Type + Genotype:Type)
dds <- DESeq(dds)
resultsNames(dds)

> coldata
                Genotype Type
C_24h_Ribo_1     Control Ribo
C_24h_Ribo_2     Control Ribo
C_24h_Ribo_3     Control Ribo
CHIK_24h_Ribo_1    Viral Ribo
CHIK_24h_Ribo_2    Viral Ribo
CHIK_24h_Ribo_3    Viral Ribo
C_24h_RNA_1      Control  RNA
C_24h_RNA_2      Control  RNA
C_24h_RNA_3      Control  RNA
CHIK_24h_RNA_1     Viral  RNA
CHIK_24h_RNA_2     Viral  RNA
CHIK_24h_RNA_3     Viral  RNA

> resultsNames(dds)
[1] "Intercept"                 "Genotype_Viral_vs_Control" "Type_RNA_vs_Ribo"          "GenotypeViral.TypeRNA"

I do not understand which contrast is the one I need, or how to access the results i'm looking for (or interpret them). Any help?

Adding to what benn pointed out, at this line:

dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~Genotype + Type + Genotype:Type)

You specified that the counts, can be a explained by Genotype (being Control or Viral) or Type (RNAseq or Ribo) or an interaction effect

What you need is this coefficient, Genotype:Type. To use it, go

res <- results(dds,name="GenotypeViral.TypeRNA")

Based on your setup, you are always comparing RNA to Ribo (see "Type_RNA_vs_Ribo" in resultsNames). If changes to the abundance of a mRNA is similar in RNAseq or Ribo, then the coefficient of the interaction term will be zero and it will be zero in the log2foldchange column in results(res). A positive fold change indicates the log2foldchange between RNA and Ribo is higher in Viral.

If you would like to compare ribo versus and rna-seq, then set RNA as the reference by:

coldata$Type = factor(coldata$Type,levels=c("RNA","Ribo"))
# run deseq again
dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~Genotype + Type + Genotype:Type)
dds <- DESeq(dds)

With this you will compare by default Ribo versus RNA, and the interaction term will tell you how much more enriched or depleted the ribo fraction is, in the Viral samples.

Would recommend checking the vignette out here:

http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#contrasts


Annotating a .vcf with centimorgan information Some programs (e.g. shapeit4) automatically annotate an INFO tag into a .vcf file which gives the cumulative genetic distance in cM between each SNP:
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##fileDate=10/11/2020 - 21:29:39
##source=G2H
##contig=<ID=9>
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
##INFO=<ID=AC,Number=1,Type=Integer,Description="Allele count">
##INFO=<ID=CM,Number=A,Type=Float,Description="Interpolated cM position">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Phased genotypes">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
##bcftools_viewVersion=1.9-207-g2299ab6+htslib-1.9-271-g6738132
##bcftools_viewCommand=view -s -1_-1 ukb_cal_chr9_v2.HumanOrigins.common_snps.HumanOriginsInds.fiftyPercentEthnicUKBInds.PHASED.bgz.vcf.gz; Date=Mon May 24 14:15:51 2021
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  -1_-1
9   257034  rs573860    A   G   .   .   AC=0;AF=0.388699;CM=0;AN=2  GT  0|0
9   270321  rs17720310  A   G   .   .   AC=0;AF=0.0820595;CM=0.007994;AN=2  GT  0|0

However, this program is really used for phasing not for .vcf annotation.
Ideally I'd like a simple command line program which takes input the .vcf and a genetic map and interpolates the genetic distance and adds a cumulative map. This could be done manually with a custom script, but quite frankly I can't be bothered to do that.
Does anyone know of a tool that exists which will do this?
I posted this to biostars a while ago and received one reply which linked to this post, which includes links to tables containing cM locations for different rs numbers.
EDIT: I am interested in doing this for human data (although the same would be the case for any diploid organism).

I think a useful tool to address your question is QCTOOL v2. Particularly, the option you are looking for is -annotate-genetic-map. This tool has been developed originally for GWAS datasets. It has multiple inputs data format.
If you could provide input data as an example, I will try working on it to provide the necessary steps to be done.
Hope this answers your question.


PLINK clump behavior on missing SNPs? I have a long list of autoimmune-associated SNPs, and I want to boil it down so that I get one SNP representing each LD block. I chose to use PLINK's --clump option for this. I'm roughly following this tutorial (but analyzing my own data).

I don't have the full data from the original GWAS's, so I am using 1000 Genomes data as the reference for LD estimation. But, the 1000 Genomes data don't have all the SNPs I need. I end up a few hundred warnings:

Warning: 'rs12345678' is missing from the main dataset, and is a top variant.

How does PLINK handle this situation? More importantly, how should I handle this situation? I don't want to omit those SNPs, but I also don't want to take them all, for the same reason that I am doing this in the first place.

Most likely, the GWASs that generated your summary statistics used other imputation panels than 1000G, like HRC. Clearly, PLINK can't estimate the LD for SNPs that aren't found in the reference, and most likely just ignores them. I think the easiest solution is to use whatever individual-level dataset you have to estimate the correlations - even if it's small or from a different population. The estimate will be biased, but you don't really need high precision for clumping.

Then again, you can skip PLINK and try searching for SNP prioritization methods based on summary statistics - it's an area of active development, and something might be available to suit your needs.


Why do NEB adapters have non-complementary sequence? I am making some infographics of different library prep types and noticed something weird about NEB's Ultra II adapters.

Molecule 1 is the desired product of the adapter ligation step. * indicates a phosphorothioate bond and - indicates a normal bond. | indicates complementary sequence through hydrogen bonding.

Molecule 2 is the product after treatment with the USER enzyme.

Molecule 3 is the result of melting molecule 2 and hybridizing the NEBNext Universal PCR Primer for Illumina.

My question is: Why does Molecule 3 have a large non-complementary region on the 3' end of the original molecule? Is that just to force Molecule 1 to have an open loop conformation so that the USER enzyme can access it better? Or maybe it is something to do with PCR kinetics?

I understand that PCR synthesizes in the 5'->3' direction, so Molecule 3 only forms with the original adapter-ligated molecules before PCR.

The original point of the non-complementary sequence is to allow formation of a Y-adapter after USER treatment, since that's what's needed for the Illumina chemistry to work. In reality the bottom strand may not be there in step 3, since USER has already digested one strand (creating the hairpin in step 2) and size selection has been performed. Regardless, it'll get removed once cluster formation is performed.


NCBI xml-format Printing nodeValues from tags with the same names I would like to retrieve all Other-source_anchor on the level of DNA binding from within a Gene Ontology XML file using php.

If I try this with my code, I only get Gene Ontology back and not the nodeValues in Other-source_anchor. Could someone tell me how to do this?

Xml file were I wwant to get the nodeValues from

<Gene-commentary_heading>GeneOntology</Gene-commentary_heading>
      <Gene-commentary_source>
        <Other-source>
          <Other-source_pre-text>Provided by</Other-source_pre-text>
          <Other-source_anchor>GOA</Other-source_anchor>
          <Other-source_url>http://www.ebi.ac.uk/GOA/</Other-source_url>
        </Other-source>
      </Gene-commentary_source>
      <Gene-commentary_comment>
        <Gene-commentary>
          <Gene-commentary_type value="comment">254</Gene-commentary_type>
          <Gene-commentary_label>Function</Gene-commentary_label>
          <Gene-commentary_comment>
            <Gene-commentary>
              <Gene-commentary_type value="comment">254</Gene-commentary_type>
              <Gene-commentary_source>
                <Other-source>
                  <Other-source_src>
                    <Dbtag>
                      <Dbtag_db>GO</Dbtag_db>
                      <Dbtag_tag>
                        <Object-id>
                          <Object-id_id>3677</Object-id_id>
                        </Object-id>
                      </Dbtag_tag>
                    </Dbtag>
                  </Other-source_src>
                  <Other-source_anchor>DNA binding</Other-source_anchor>

code:

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

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

use a XSLT processor (http://php.net/manual/en/class.xsltprocessor.php) with the following stylesheet.

<?xml version='1.0' encoding="UTF-8"?>
<xsl:stylesheet xmlns:xsl='http://www.w3.org/1999/XSL/Transform' version='1.0'>
    <xsl:output method="html"/>

    <xsl:template match="/">

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

    </xsl:template>

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

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

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

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

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

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

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

    </xsl:template>

</xsl:stylesheet>


Is it possible to check if a patient has the HLA-B27 antigen from his exome stored in a VCF file Is it possible to check if patient has the HLAB27 antigen by comparing his exome stored in VCF with all known HLAB27 variants?

For instance, rs4349859 is a known HLAB27 SNP.

If I can find rs4349859 in the patient's VCF file, could I say he is HLAB27+ ?

Some papers claim that you can do HLA typing from exome sequencing data. I have not validated this myself. I also guess it depends on capture. Nonetheless, I don't think you can type HLA with VCF generated with a standard pipeline, which is not specifically designed for HLA. Current variant callers won't work on HLA genes due to their high diversity. Also in general, there are thousands of HLA-B alleles. In most cases, you can't find the allele by looking at only one SNP.

No, you can't check if the patient has HLA-B27 only from VCF. There might be a hope if you get raw reads and run a proper HLA typing tool designed for such data.


How to Installing & Using RDP Classifier through bioconda? I am interested in using RDP classifier and RDPTools at the command line.  In following the docs from RDP, I ran into this issue during the build, after running make:

BUILD FAILED
/home/denir/bin/RDPTools/Clustering/nbproject/build-impl.xml:896: The following error occurred while executing this line:
/home/denir/bin/RDPTools/Clustering/nbproject/build-impl.xml:1470: The following error occurred while executing this line:
/home/denir/bin/RDPTools/AlignmentTools/nbproject/build-impl.xml:896: The following error occurred while executing this line:
/home/denir/bin/RDPTools/AlignmentTools/nbproject/build-impl.xml:1414: The following error occurred while executing this line:
/home/denir/bin/RDPTools/ReadSeq/nbproject/build-impl.xml:920: The following error occurred while executing this line:
/home/denir/bin/RDPTools/ReadSeq/nbproject/build-impl.xml:300: Compile failed; see the compiler error output for details.

Total time: 1 second
Makefile:15: recipe for target 'Clustering/dist/Clustering.jar' failed
make: *** [Clustering/dist/Clustering.jar] Error 1

As one user suggested, I have tried downgrading from ant 1.10 to ant 1.9, this didn't work for me.

I would rather not downgrade to JDK 7 as another user suggested if possible so I am looking for alternatives. One user was able to install rdp classifier with bioconda. I was able to do that OK. However, I'm not sure where to go from here? The tutorial provided with commands like:

java -Xmx2g -jar ~/RDPTools/Clustering.jar derep -m '#=GC_RF' -o derep.fa all_seqs.ids all_seqs.samples alignment/*.fasta

Depend on the build created using the install method recommended by RDP.

Does anyone else use rdp classifier installed via conda?  Is there a quickstart guide out there that I could refer to?  Thanks so much for your help!

The benefit of using the conda package is that you no longer need to know or care exactly where the jar files are (that many java-based packages require this is TERRIBLE design). The example java -jar ... command in your post becomes:

Clustering derep -m '#=GC_RF' -o derep.fa all_seqs.ids all_seqs.samples alignment/*.fasta

in the bioconda version. In bioconda packages, there are wrapper (usually python) scripts meant to provide somewhat friendlier interfaces to java programs that otherwise lack executables. In this case, RDP Classfier has the following executables in bioconda:

AbundanceStats
AlignmentTools
Clustering
FrameBot
KmerFilter
ProbeMatch
ReadSeq
SeqFilters
SequenceMatch
hmmgs
classifier

If you really want to know, the .jar file is under $PREFIX/share/rdptools-$PKG_VERSION-$PKG_BUILDNUM/, where $PREFIX is where ever the environment is and the other two variables indicate the package version and build number (at the moment, rdptools-2.0.2-1 would be the likely folder name).


Installing RnBeads via bioconductor - .onLoad failed in loadNamespace(), call: NULL I am trying to install RnBeads from bioconductor but the installation of its dependancy TxDb.Hsapiens.UCSC.hg19.knownGene fails. The error is:

biocLite("RnBeads")
BioC_mirror: https://bioconductor.org
Using Bioconductor 3.6 (BiocInstaller 1.28.0), R 3.4.3 (2017-11-30).
Installing package(s) ‘RnBeads’
also installing the dependencies ‘FDb.InfiniumMethylation.hg19’, ‘methylumi’

trying URL 'https://bioconductor.org/packages/3.6/data/annotation/src/contrib/FDb.InfiniumMethylation.hg19_2.2.0.tar.gz'
Content type 'application/x-gzip' length 38996322 bytes (37.2 MB)
==================================================
downloaded 37.2 MB

trying URL 'https://bioconductor.org/packages/3.6/bioc/src/contrib/methylumi_2.24.1.tar.gz'
Content type 'application/x-gzip' length 6178909 bytes (5.9 MB)
==================================================
downloaded 5.9 MB

trying URL 'https://bioconductor.org/packages/3.6/bioc/src/contrib/RnBeads_1.10.2.tar.gz'
Content type 'application/x-gzip' length 6437974 bytes (6.1 MB)
==================================================
downloaded 6.1 MB

* installing *source* package ‘FDb.InfiniumMethylation.hg19’ ...
** R
** data
** demo
** inst
** preparing package for lazy loading
Error: package or namespace load failed for ‘TxDb.Hsapiens.UCSC.hg19.knownGene’:
 .onLoad failed in loadNamespace() for 'TxDb.Hsapiens.UCSC.hg19.knownGene', details:
  call: NULL
  error: evaluation nested too deeply: infinite recursion / options(expressions=)?
Error : package ‘TxDb.Hsapiens.UCSC.hg19.knownGene’ could not be loaded
ERROR: lazy loading failed for package ‘FDb.InfiniumMethylation.hg19’

The downloaded source packages are in
    ‘/tmp/RtmpY8Yf4l/downloaded_packages’
installation path not writeable, unable to update packages: doMC, doSNOW,
  foreach, iterators, snow, boot, class, cluster, codetools, foreign,
  KernSmooth, lattice, MASS, Matrix, mgcv, nlme, nnet, rpart, spatial, survival

As suggested by the error message, I tried setting options(expressions=500000), which unfortunately did not change the error message. An important note may be that I am under network mount and don't have administrative access rights.

Edit

I don't have a .Rprofile file.

Note about installation path not writeable:
Typing .libPaths() reveals that my home R folder is at the first position (and it is writable). The other paths ("/usr/local/lib/R/site-library"                          , "/usr/lib/R/site-library" and "/usr/lib/R/library") are not writable without sudo.

I also tried to delete my home R folder and to reinstall it. The reinstallation of all packages worked for all packages except for the three listed above (same error message). I also tried biocLite("TxDb.Hsapiens.UCSC.hg19.knownGene", lib.loc=.libPaths()[1], lib=.libPaths()[1]) to specifically address my home R directory, but the same error is showing up (but without the installation path not writeable part).

Here is my sessionInfo():

R version 3.4.3 (2017-11-30)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.4 LTS

Matrix products: default
BLAS: /usr/lib/openblas-base/libopenblas.so.0
LAPACK: /usr/lib/lapack/liblapack.so.3.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
 [9] LC_ADDRESS=C               LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

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

other attached packages:
[1] BiocInstaller_1.28.0

loaded via a namespace (and not attached):
[1] compiler_3.4.3 tools_3.4.3

Thank you very much!

I have reported the error on support.bioconductor.org but we were not able to solve the problem. It seems that some dependencies are making problems, but as reinstalling my home R folder did not resolve the issue, I have come up with the following workaround:
I downloaded the latest R version from my favorite CRAN mirror, opened an R session with pathToNewR/bin/R from my terminal and installed RnBeads there. This works fine.
Note: I had to make sure that pathToNewR/library is the only entry in .libPaths(), as per default, my home R folder was also contained in .libPaths().
As a consequence, I am independent from the (maybe broken) R version and packages on my server and can finally work with RnBeads :)


Filtering bam file with bedtools intersect using strand specificity I have a small sample of records from a bam file. I obtained them by intersecting with a bed file not enforcing strand-specificity. So I get them all back upon performing the same intersection again:

$ bedtools intersect -a sample.bam -b features.bed -wa -u -f 1.0 | samtools view
MN00481:17:000H2FH5N:1:11104:12480:15984    16  III 4499483 42  28M *   0   0   CAGGAACACGCACATCGTCGATTTTCAG    FFFFFAFFFFFFFFFFFFFFFFFFFFFF    AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:28 YT:Z:UU
MN00481:17:000H2FH5N:1:11103:9925:11522 16  IV  2388448 42  30M *   0   0   CACTGGAACTGGTCCGATCTTCAAGCCGGT  FFFFFFFFFFA/FFFFAFFFFFF/FFAFAF  AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:30 YT:Z:UU
MN00481:17:000H2FH5N:1:22109:11939:5243 16  IV  4749336 42  29M *   0   0   ATGTGGATGCGGTCAGCTTGGGCACGATG   FFFFFFFFFFFF/FFFF6FFFFFFFF/FF   AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:29 YT:Z:UU
MN00481:17:000H2FH5N:1:22104:17385:19431    16  V   4443064 42  29M *   0   0   CAAAATGTGAACAGGAACGTTTTCATGAG   F/FFFFFFFFFFFFFFFFFFFFFFFFFFF   AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:29 YT:Z:UU
MN00481:17:000H2FH5N:1:23105:22058:12968    0   X   981524  42  29M *   0   0   CCATCGATTCTCATTTTAATGCTGGTGCT   FFFFFFFFFFFFFFFFFFFFFFFFFFFFF   AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:29 YT:Z:UU
MN00481:17:000H2FH5N:1:12104:7867:7500  16  X   4214381 42  29M *   0   0   GCTGCAGGAATGACGGCTCTAGCAAACGT   FFFFFFFFFFFFFFFFFFFFFFFFFFFFF   AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:29 YT:Z:UU

Now I would like to enforce strand-specificity, but in either case, none of the records pass the filter:

Same strand (-s):

$ bedtools intersect -a sample.bam -b features.bed -wa -u -f 1.0 -s | samtools view
$ # nothing here

Opposite strand (-S):

$ bedtools intersect -a sample.bam -b features.bed -wa -u -f 1.0 -S | samtools view
$ # nothing here

What is the correct way of doing this?

Edit

Bedtools version:

$ bedtools --version
bedtools v2.27.1

My bed file has strand information:

$ tail features.bed
X   17617826    17619458    "WBGene00015867"    +   .
X   17622777    17640743    "WBGene00019211"    -   .
X   17650177    17651594    "WBGene00044771"    +   .
X   17662457    17665357    "WBGene00003744"    +   .
X   17671028    17680790    "WBGene00020732"    +   .
X   17680910    17682294    "WBGene00044772"    -   .
X   17695231    17701744    "WBGene00017625"    +   .
X   17710485    17711674    "WBGene00019189"    -   .
X   17714455    17714832    "WBGene00007069"    +   .
X   17714967    17718717    "WBGene00007068"    +   .

Your bed file switch the score and strand columns. Strand should be column 6, score column 5. For example take this bed file with the 5th and 6th columns switched - zero results for strand specific intersection b/c strand is actually '.'

$ cat wrong.bed
chr1    30000000        31000000        noname  +       .
$ bedtools intersect -a mm10_500k_pe.bam -b wrong.bed -wa -u -f 1.0 \
    | samtools view | wc -l
166
$ bedtools intersect -a mm10_500k_pe.bam -b wrong.bed -wa -u -f 1.0 -s \
   | samtools view | wc -l
0

After fixing the bed file:

$ cat right.bed
chr1    30000000        31000000        noname  .       +
$ bedtools intersect -a mm10_500k_pe.bam -b right.bed -wa -u -f 1.0 \
      | samtools view | wc -l
166
$ bedtools intersect -a mm10_500k_pe.bam -b right.bed -wa -u -f 1.0 -s \
      | samtools view | wc -l
83


How to convert FASTA to BED I have a FASTA file:

> Sequence_1
GCAATGCAAGGAAGTGATGGCGGAAATAGCGTTAGATGTATGTGTAGCGGTCCC...
> Sequence_2
GCAATGCAAGGAAGTGATGGCGGAAATAGCGTTAGATGTATGTGTAGCGGTCCC....
....

I want to generate a BED file for each sequence like:

Sequence_1 0 1500
Sequence_2 0 1700

The BED regions will simply be the size of the sequences.

Q: I did that before with a one-line command. I don't remember what that is, it was on Biostars. I can't find the post now. What's the simplest way to do the conversion?

You can do this easily with bioawk, which is a version of awk with added features facilitating bioinformatics:

bioawk -c fastx '{print $name"\t0\t"length($seq)}' test.fa

-c fastx tells the program that the data should be parsed as fasta or fastq format. This makes the $name and $seq variables available in the awk commands.


How do I create a VCF file of all known pathogenic mutations in a gene of interest? I would like to create a list (probably .vcf format would be good) of all known pathogenic missense mutations in a human gene of interest and then add other variants that could lead to the same pathogenic amino acid substitutions. I intend to search targeted sequencing data of my gene of interest for variants that would cause these particular substitions.
I'm very new to bioinformatics so I'm not so familiar with the appropriate databases and tools to use, but this is a rough idea of how it could be done:

extract all mutations listed in Ensembl within my gene of interest that are labelled with HGMD_MUTATION  (is there an option to download a vcf file with all HGMD mutations within my gene? I've only managed to do so for individual mutations.)
annotate list with reference and mutant amino acids for every mutation in the list
add additional mutations (which may not be present in HGMD) which could lead to the same substitution. E.g. I'm interested in the substitution Proline (reference codon: CCG) to Leucine (TTR or CTN) at AA position 102, so I would like to list all hypothetical mutations that would create a leucine codon in this position, not only the one given in Ensembl (which is CCG -> CTG).

Could you point me towards the best tools to use for these purposes? Unfortunately I haven't got a HGMD subscription. Many thanks!

I'm not sure how to generate the additional mutations, but I would say that HGMD is not the way to find all the pathogenic variants. I would probably filter this table by either Clinical significance or Evidence->Phenotype association.


Run gffread in multi-thread mode Is there any possibility to run gffread in multi-thread mode? The answer seems to be 'no' from the manual (or gffread -h), as no multi-thread option is mentioned.

I'm mostly using this utility to extract transcript sequences (FASTA) from annotation files (GTF). In single-thread mode, runtime is acceptable in most cases, but would like to speed up a bit in some cases (especially for big genomes).

Note that I already tried, as mentioned in the manual, to index the genome first (which, again, cuts the runtime in most but not all cases):

  Note that the retrieval of the transcript sequences this way is going to be much faster if a fasta index file (genome.fa.fai in this example) is found in the same directory with the genomic fasta file.

  Such an index file can be created with the samtools utility prior to running gffread, like this: samtools faidx genome.fa

  Then in subsequent runs using the -g option gffread will find this
  FASTA index and use it to speed up the extraction of transcript
  sequences.

I'm also open to different alternatives than gffread to convert GTF -> FASTA.

Since there doesn't seem to be an easy way to run this in parallel, you could instead break the job into sections. For example, separate each chromosome into its own gtf file, extract the sequences using that file and then cat them all together.

The commands would be something like this (using the human hg38 genome and the gff GENCODE annotations downloaded from here):

Split the file into one file per chromosome (or whatever else you have as the first field of your gff/gtf file):

tmp=$(mktemp -d); ## create a temp directory
awk -vtmp="$tmp" '($1!~/^#/){print > tmp"/"$1".gff"}' gencode.v27.annotation.gtf

Extract the sequences described in each file, running each command in the background (&):

for f in $tmp/*; do
    gffread -w "$f.fa" -g hg38.fa "$f" &
done

Concatenate them into one file:

cat "$tmp/*fa" > all.fa

You can combine all this into a single command, using wait to make sure you don't concatenate until all of the sub-commands have finished:

tmp=$(mktemp -d);
awk -vtmp="$tmp" '($1!~/^#/){print > tmp"/"$1".gff"}' gencode.v27.annotation.gtf &&
    for f in "$tmp"/*; do
        gffread-0.9.12.Linux_x86_64/gffread -w $f.fa -g hg38.fa "$f" &
    done
    wait
    cat "$tmp"/*fa > all.fa

Or, just save it as a script called pgffread.sh (or whatever):

#!/bin/sh
gtf="$1"
genome="$2"

tmp=$(mktemp -d);
awk -vtmp="$tmp" '($1!~/^#/){print > tmp"/"$1".gff"}' "$gtf" &&
    for f in "$tmp"/*; do
      gffread -w "$f".fa -g "$genome" "$f" &
    done
wait
cat "$tmp"/*fa
rm -rf "$tmp"

And then run it giving the genome and annotation files as parameters:

pgffread.sh gencode.v27.annotation.gtf hg38.fa > all.fa

On the system I ran it on, which has 48 cores, the time difference was significant (although it's so fast anyway, I really don't know if it is worth the effort unless you do this very often):

$ time gffread -w out.fa -g hg38.fa gencode.v27.annotation.gtf
real    1m52.072s
user    1m3.160s
sys     0m47.257s

$ time ./pgffread.sh gencode.v27.annotation.gtf hg38.fa > all.fa
real    0m17.533s
user    1m5.228s
sys     0m43.039s


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

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

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


Cellranger aggr and Seurat merge difference? I have multiple libraries of 10x Chromium single-cell RNA-seq data, which I'd like to combine. One option is using cellranger aggr which by default does a depth normalization:

mapped: (default) Subsample reads from higher-depth libraries until they all have an equal number of confidently mapped reads per cell. raw: Subsample reads from higher-depth libraries until they all have an equal number of total (i.e. raw, mapping-independent) reads per cell.
none: Do not normalize at all.

https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/aggregate

The new version of Seurat (2.0) provides another way via the MergeSeurat() (or AddSamples()) functions.

http://satijalab.org/seurat/merge_vignette.html

These have the arguments:

do.logNormalize whether to normalize the expression data per cell and transform to log space.
total.expr      scale factor in the log normalization
do.scale        In object@scale.data, perform row-scaling (gene-based z-score)
do.center       In object@scale.data, perform row-centering (gene-based centering)

I'm wondering how do these methods compare? What are the pros/cons etc? What do we need to consider when we combine samples?

Cell Ranger aggregate subsamples reads (unless you select none), so you will end up with less total reads in samples that have more initially. The output is still raw counts, but you will have more or less per cell.

Seurat just merges the raw counts matrices and normalizes those.


Is there an efficient way to check an input BAM in R? I'm writing a function in R for an R package which takes as input a BAM.

my_func = function(my_bam){
    if (!file.exists(my_bam)){
        stop("Input a valid BAM file.")
    }

    ## do awesome stuff with BAM

}

Trying to write robust software, I would like to do more of a check whether this BAM is actually a BAM. Above, I'm just using file.exists().

I could also do something simple like "check whether the file has the extension *.bam" with

grepl('*.bam', my_bam)  ## TRUE for file.bam, FALSE for file.txt

However, is there a quick way to check that the input is actually a BAM? I would prefer to throw an error here, if the user inputs a *txt or *bedpe.

I don't think it's necessary to validate the BAM per se, whereby I believe most users rely upon Picard:

https://broadinstitute.github.io/picard/command-line-overview.html#ValidateSamFile

My first thought was Rsamtools for this operation:

https://bioconductor.org/packages/release/bioc/html/Rsamtools.html

However, it doesn't appear that there is a function here which checks the user inputs a BAM, e.g. BamFile() doesn't appear to throw an error

> library(Rsamtools)
> fake_bam = 'nonsense_string.bam'   ## a string I made up, not a real BAM
> BamFile(fake_bam)
class: BamFile
path: nonsense_string.bam
index: NA
isOpen: FALSE
yieldSize: NA
obeyQname: FALSE
asMates: FALSE
qnamePrefixEnd: NA
qnameSuffixStart: NA

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

You won’t know whether it’s a valid BAM file until you attempt reading it in full; even Devon’s method (which I recommend) only does a superficial check — i.e. it checks whether (parts of) the BAM header can be read correctly. However, this doesn’t guarantee that the file isn’t truncated or contains invalid data.

This is of course a more general problem, but it’s perhaps a bit stranger for BAMs (and other biological formats), as it’s relatively frequent to end up with technically invalid files that were produced by existing software. The only recourse is actually trying to use the file, and to handle errors appropriately.

For a very quick check™ you could also pass the file through gunzip and test whether the “magic bytes” at the beginning of the file equal 'BAM\1':

bgzf_stream = gzfile(filename, 'r')
magic = readChar(bgzf_stream, 4)
stopifnot(identical(magic, 'BAM\1'))


How do I rewrite a read group using pysam? I am trying to rewrite a SAM/BAM file with altered read group entries using pysam. In this simplified version, I want to take a BAM, and rewrite the SM tags in all read groups to the same string and copy the alignments with new header into a new BAM. The contents of the read group dictionaries I am changing appear to have no effect on the header that prints or is written to file output. How can I alter the read groups?

import argparse
import pysam

def rewrite_read_groups(input_filename, output_filename, sample):
    try:
        bam_input = pysam.AlignmentFile(input_filename, 'rb')
        # replace header fields
        output_header = bam_input.header
        read_groups = output_header['RG']
        print(read_groups)
        for read_group in read_groups:
            read_group['SM'] = sample
            print(read_group)
        output_header.to_dict()
        print(output_header)
    finally:
        bam_input.close()

    print(output_header)

    try:
        bam_input = pysam.AlignmentFile(input_filename, 'rb')
        bam_output = pysam.AlignmentFile(output_filename, 'wb', header=output_header)

        # copy alignments
        for read in bam_input.fetch(until_eof=True):
            bam_output.write(read)

    finally:
        bam_input.close()
        bam_output.close()

if __name__ == "__main__":
    parser = argparse.ArgumentParser(description="Rewrite bam read groups")
    parser.add_argument('-s', "--sample", help="Set sample for all read groups")

    parser.add_argument('-i', "--input_file", help="Input bam", required=True)
    parser.add_argument('-o', "--output_file", help="Output bam", required=True)

    args = parser.parse_args()

    rewrite_read_groups(args.input_file, args.output_file, args.sample)

If the goal is to assign all alignments to a single read group (which this code seems to do), then Picard AddOrReplaceReadGroups might help: https://broadinstitute.github.io/picard/command-line-overview.html#AddOrReplaceReadGroups.

Modifying headers using PySAM might be tricky – the internal representation of headers has recently been changed and seems not to be stable across PySAM versions.


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

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

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

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


How to efficiently extract coordinates from a `pybedtools.Interval`? I'm using pybedtools to process data, and at some point, I loop over a BedTool object and use unpacking to extract coordinates from the Interval:

for interval in my_bedtool:
    (chrom, start, end, _, _, strand) = interval

Suppose I want to compute the interval length, as follows:

    length = end - start

This results in a TypeError: unsupported operand type(s) for -: 'str' and 'str'.

It happens that the unpacked end and start variables are strings.
Converting them to ints is probably suboptimal given that Intervals can be created by providing ints for the coordinates:

from pybedtool import Interval
interval = Interval("X", 0, 100, strand="+")

The start and end attributes of an Interval are ints.

I could do:

    (chrom, start, end, strand) = (
        interval.chrom, interval.start, interval.end, interval.strand)

But I think this looks a bit ugly.

Is there a better approach?

The reason why unpacking an Interval results in strings instead of ints might be related to how indexing by position works.

Even without considering the string/int issue, it seems slightly more efficient to access the interval's attributes rather than unpacking them directly:

In [311]: %timeit [(chrom, start, end, strand) for (chrom, start, end, _, _, strand) in my_bedtool]
306 ms ± 26.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [312]: %timeit [(interval.chrom, interval.start, interval.end, interval.strand) for interval in my_bedtool]
238 ms ± 24.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

To resolve the ugliness issue, attribute access can be made more compact through the use of operator.attrgetter:

In [313]: from operator import attrgetter

In [314]: get_interval_attrs = attrgetter("chrom", "start", "stop", "strand")

In [315]: %timeit [get_interval_attrs(interval) for interval in my_bedtool]
239 ms ± 26.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)


merging two/or more bed file into one bed file I am trying to merge two bed files (more in future) to one. my bed files are something like :
.

I need to merge them in a way to have the shared chromosome location.

Is there a way to do that ?

In my opinion, the easiest way to merge bed files is to use bedtools merge.


Why does the Genbank phylogenetic tree of 4000 sars-cov-2 sequences only display 200 nodes? I manually selected all 4000 sequences on the 20 pages of the SARS-CoV-2 sequences on genbank here.
When I click on Build phylogenetic tree, it only has 200 nodes.
Is this expected? I was expecting to see 4000 nodes.
Is it simply a displaying issue where only 200 sequences per page are revealed while browsing the list of sequences? Is there a way to get the tree for all 4000 sequences?

Okay, I'll provide a formal answer, Genbank phylogeny is to get a gist of the diversity. It isn't a formal phylogenetic analysis. You can down the tree as a pdf, I don't think you can download the treefile (in fact I'm pretty certain you can't). Genbank will perform a reasonable distance method and from recollection it will "collapse clades", so if you see little triangles on the tree it represents a large number of taxa represented as a triangle.

I was surprised to learn there are 4000 SARS-CoV-2 sequences.

To make a tree you need to download the data, align (muscle) and send through a formal phylogeny program. If you think Genbank is okay, then try MegaX its very point and clicky. Generally, maximum likelihood phylogeny has returned to vogue for at least 1 decade now (previously Bayesian). However, the Beast package is a Bayesian analysis which is still going strong and the main tool now used in molecular dating.

In summary, to perform a tree analysis beyond Genbank

use the download option (there are two variants - the entire genome or the cut locus based on the query sequence)
align using muscle or clustalo (muscle is really quick(
perform a tree analysis using MEGA X, phylml, raxml or Beast depending on what you objectives are


Which software is used to run Molecular Dynamics simulation of a small RNA hairpin? I am currently working in a project that needs MD simulation of a small RNA hairpin. I am currently trying it on GROMACS but it seems to be protein based, hence am facing several issues.
Is there any software that is designed for the same?

Google found this review that might be of interest: https://pubs.acs.org/doi/10.1021/acs.chemrev.7b00427

Here is a link to a free full-text version:
https://www.ncbi.nlm.nih.gov/pmc/articles/pmid/29297679/

The authors "here provide a comprehensive overview of the fast-developing field of MD simulations of RNA molecules. We begin with an in-depth, evaluatory coverage of the most fundamental methodological challenges that set the basis for the future development of the field, in particular, the current developments and inherent physical limitations of the atomistic force fields and the recent advances in a broad spectrum of enhanced sampling methods. We also survey the closely related field of coarse-grained modeling of RNA systems. After dealing with the methodological aspects, we provide an exhaustive overview of the available RNA simulation literature, ranging from studies of the smallest RNA oligonucleotides to investigations of the entire ribosome. Our review encompasses tetranucleotides, tetraloops, a number of small RNA motifs, A-helix RNA, kissing-loop complexes, the TAR RNA element, the decoding center and other important regions of the ribosome, as well as assorted others systems. Extended sections are devoted to RNA–ion interactions, ribozymes, riboswitches, and protein/RNA complexes. Our overview is written for as broad of an audience as possible, aiming to provide a much-needed interdisciplinary bridge between computation and experiment, together with a perspective on the future of the field."

Hope this helps.


How can HISAT2/StringTie report decimal coverage values I have performed RNA-seq analysis using HISAT2 & StringTie workflow suggested in: Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown.

Some of the transcripts/exons have decimal coverage values (eg., 1.1 or 2.59).

My question is: How can these tools report decimal coverage? If only half of the read overlaps exon will this exon have 0.5 coverage?

I suspect this is average coverage across the transcript or exon.

Formula for average coverage:

Where:

N: number of reads mapping/aligning to your transcript/exon
L: average read length
T: Transcript/exon length


How to calculate logCPM across all samples? Using edgeR for differential analysis between Tumor and Normal gave me differential expressed genes with logFC, logCPM, PValue and FDR.

From the details of glmTreat function I see that logCPM is average log2-counts per million, the average taken over all libraries.

The output table I got from glmTreat:

Geneid     logFC  unshrunk.logFC    logCPM  PValue       FDR
Gene1   4.985298235 5.573931983 1.236660486 1.55E-54    1.87E-50
Gene2   4.613438634 5.52126484  0.920343233 2.59E-53    1.56E-49
Gene3   5.250601296 5.356432653 2.751294034 2.48E-50    9.95E-47
Gene4   4.943049159 5.379741165 1.361393757 2.06E-45    6.19E-42
Gene5   6.100121754 6.117580436 5.774385315 2.95E-43    7.11E-40
Gene6   4.722697891 5.320461275 1.120685402 3.05E-42    6.11E-39
Gene7   5.246129853 5.497012001 1.902992053 2.78E-40    4.78E-37
Gene8   3.878773277 4.956636276 0.776208741 8.66E-39    1.30E-35
Gene9   4.441752496 4.94930499  1.132652682 4.37E-38    5.83E-35

To calculate logCPM given in the above table manually I did like this:

logCPM <- cpm(y, prior.count=2,  log=TRUE)

I got the logCPM values for all the genes. Then I calculated Average across all the samples for each gene. I got like below:

Geneid  Average
Gene1   0.686560246
Gene2   0.617115826
Gene3   1.075975225
Gene4   0.692050878
Gene5   1.277556065
Gene6   0.638358189
Gene7   0.689323163
Gene8   0.60700396
Gene9   0.662115092

Why there is so difference? Where I'm doing wrong?

If you want the values by group then subset y to contain the samples of interest and then feed that to aveLogCPM().


Calling GEOquerry error while loading the library I started getting a strange error when trying to load GEOquerry with library(GEOquerry).
Error: package or namespace load failed for ‘GEOquery’ in loadNamespace(i, c(lib.loc, .libPaths()), versionCheck = vI[[i]]): namespace ‘rlang’ 0.4.9 is already loaded, but >= 0.4.10 is required

When I run > devtools::install_github('r-lib/rlang') I get ...
Error in loadNamespace(i, c(lib.loc, .libPaths()), versionCheck = vI[[i]]) :
  namespace ‘rlang’ 0.4.9 is already loaded, but >= 0.4.10 is required

Has anyone dealt with this before?  Restarting my R session has no effect.

rlang 0.4.10 is the current version at CRAN: https://cran.r-project.org/web/packages/rlang/index.html
Try to install it via install.packages, not via GitHub directly.


scRNASeq expression matrix with decimal values I am trying to replicate some results of a scRNASeq experiment and, when I looked at the data provided by the author, I noticed that some of the counts in the expression matrix are represented as decimals.

Would this be because the authors uploaded a normalised count matrix or for other reasons? Also, please note that the authors used the 10X genomics cellranger pipeline for the processing of raw data.

I have provided a brief snippet of the data frame here:

structure(list(P1TLH_AAACCTGAGCAGCCTC_1 = c(0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), P1TLH_AAACCTGTCCTCATTA_1 = c(0,
0, 0.314759545319035, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0.314759545319035, 0, 0), P1TLH_AAACCTGTCTAAGCCA_1 = c(0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1.18024827050334, 0, 0, 0, 0, 0, 1.18024827050334,
0, 0, 0), P1TLH_AAACGGGAGTAGGCCA_1 = c(0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0.644905372424165, 0, 0, 0), P1TLH_AAACGGGGTTCGGGCT_1 = c(0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.780522072779651,
0, 0, 0), P1TLH_AAAGCAACAGTAAGAT_1 = c(0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)), row.names = c("RP11-34P13.7",
"FO538757.2", "AP006222.2", "RP4-669L17.10", "RP5-857K21.4",
"RP11-206L10.9", "LINC00115", "FAM41C", "RP11-54O7.1", "RP11-54O7.3",
"SAMD11", "NOC2L", "KLHL17", "PLEKHN1", "RP11-54O7.17", "HES4",
"ISG15", "AGRN", "C1orf159", "RP11-465B22.8"), class = "data.frame")

Also, I am quite confused of the structure they used - as I don't have much experience with single cell datasets. Just to give some further context:

The colnames of the "full" dataframe are coded from P1TLH_sequence_1 to P5TLH_sequence_1. There are 5 patients in the study.

I'm assuming that each of the colnames represents a cell for a particular patient while the rows are the genes.

If this assumption is correct, how would I do the conversion from colnames in the dataframe to the "cell names"? And lastly, why are there decimals in this expression matrix?

Thank you all for your time,

It looks like you have normalized counts. My guess is that they are log-normalized with some scale factor, probably 10000.

With the whole matrix you should be able to work out the scale factor for each cell. See https://www.nxn.se/valent/2018/10/25/unscaling-scaled-counts-in-scrna-seq-data

They typical format for scRNA-seq data is to have cells in the columns and genes in the rows. So here your rownames are the gene names and colnames are the cell names.


How to superimpose barplots in karyotype graphic using the chromPlot Bioconductor package? I have a data frame containing 4 columns:

Chrom: chromosome name of the feature
Start: starting base of the feature
End: end base of the feature
Name: name of the feature, e.g. 5S ribosomal RNA

The head of the data frame is shown below:

      Chrom    Start      End                        Name
    1     1 15390242 15390720 16S ribosomal RNA (partial)
    2     1 15390900 15391794 16S ribosomal RNA (partial)
    3     1 15392798 15396046           23S ribosomal RNA
    4     1 15397439 15397917 16S ribosomal RNA (partial)
    5     1 15398097 15398991 16S ribosomal RNA (partial)
    6     1 15399995 15403243           23S ribosomal RNA

I want to plot these features on the karyotype of C. elegans. I am using the chromPlot package from Bioconductor.

I am using the following code:

Subsetting the original data frame to get a data frame for each feature type: 5S, 16S and 23S rRNAs

annot1 <- dplyr::filter(cel, Name == '5S ribosomal RNA')
annot2 <- dplyr::filter(cel, Name == '23S ribosomal RNA')
annot3 <- dplyr::filter(cel, Name == '16S ribosomal RNA (partial)')

and then plot:

chromPlot(gaps = cel, figCols = 2, annot1 = annot1, annot2 = annot2)

This is OK because the two features are far away in the genome. However, if I try:

chromPlot(gaps = cel_gap, figCols = 2, annot1 = annot1, annot2 = annot2,
          annot3 = annot3)

one of the features masks the other because they are close to each other in the genome. How can I make the green bar stack on top of the yellow one instead of covering it?

chromPlot(gaps = cel_gap, figCols = 2, annot1 = annot1, annot2 = annot3, annot3 = annot2,
          chrSide=c(1, 1, -1, 1, -1, 1, -1, 1))

The chrSide argument decides on which side the histogram goes, so it can be used to solve my problem. The bars are not stacked but the result is the same: more accurate visualisation.

Thanks to Llopis for the guidance.