How can I calculate gene_length for RPKM calculation from counts data? I have read counts data and I want to convert them into RPKM values. For this conversion I need the gene length.

Does the gene length need to be calculated based on the sum of coding exonic lengths? Or are there any different ways for that?

I know that gene length can be taken from the Gencode GTF v19 file. Could you please tell me how that Gene_length is calculated?

The Data I'm having is RNA-Seq data. I don't have any idea whether I need to include UTR's in this calculation or only exons?

In Github I have seen RPKM calculation from Counts data with the Gene_length from Gencode GTF file. Do you think this is the right way of calculation?

And why RPKM is - Its not for differential analysis. For TNBC subtyping they use microarray data. I would like to give a try with RNA-Seq data. So for this I'm trying out different and the right way.

Here you can find some example R code to compute the gene length given a GTF file (it computes GC content too, which you don't need). This uses one of a number of ways of computing gene length, in this case the length of the "union gene model". In this method, the non-duplicated exons for each gene are simply summed up ("non-duplicated" in that no genomic base is double counted). This is a very simple way of getting a gene length.

There are alternative methods that you should be aware of, among which are:

Median transcript length: That is, the exonic lengths in each transcript are summed and the median across transcripts is used. This is probably a little more valid than the code that I linked to.
Per-sample effective gene lengths: the optimal method, though it requires using something like RSEM, which will give you an effective gene length.
Using the length of the "major isoform" in your tissue of interest. This isn't as good as method 2, but is more accurate than all of the others.

At the end of the day, you're just coming up with a scale factor for each gene, so unless you intend to compare values across genes (this is problematic to begin with) then it's questionable if using some of the more correct but also more time-involved methods are really getting you anything.

Edit: Note that if you want to plug these values into some sort of subtyping tool (TNBC in your case), you should first start with some samples for which you know the subtype. Then you can at least see if you're getting reasonable results. After that, do read up on how the method works and see if there's anything about RNAseq that makes it incompatible.


Fast processing of fastq data I am trying to write python script for customized filtering for fastq file (size >3 GB). My proposed script is as follows:

def filtering(read):
    time.sleep(0.1)

    if len(read) >= 15 and \
       len(read) <= 30 and \
       np.mean(read.qual) >= 30 and \
       sum(i < 20 for i in read.qual) <=2:
           flag1 = True
           return  flag1

pool = multiprocessing.Pool(5)
pool.ncpus = 4
for file in os.listdir('/moRNA_data/trimmed_fastq'):
    if file.endswith('.fq'):
        fastq_file = HTSeq.FastqReader( os.path.join('/moRNA_data/trimmed_fastq',file))
        name = str(file.split('.')[0]) + '.fastq'
        file_out = open(os.path.join('/moRNA_data/filtered_fastq',name), "w" )
        for read in tqdm(fastq_file):
            flag1 = pool.apply_async(filtering,read)
            if flag1:
                read.write_to_fastq_file( file_out )
pool.close()
pool.join()

The problem with this script is that it works well with small file like upto 100-200 MB but when it goes to directory having multiple sample each having size more than 3
GB, it works well for first file but after that it hangs so severely that I had to restart the server [ I can't even kill the process]. I need suggestions to improve this script to make fastq file processing faster. Platform details is as follows:

OS : Ubuntu-18.02
Architecture:        x86_64
CPU op-mode(s):      32-bit, 64-bit
Byte Order:          Little Endian
CPU(s):              12
On-line CPU(s) list: 0-11
Thread(s) per core:  2
Core(s) per socket:  6
Socket(s):           1
NUMA node(s):        1
Vendor ID:           GenuineIntel
CPU family:          6
Model:               158
Model name:          Intel(R) Core(TM) i7-8700
CPU @ 3.20GHz
Stepping:            10
CPU MHz:             800.560
CPU max MHz:         4600.0000
CPU min MHz:         800.0000
BogoMIPS:            6384.00
Virtualization:      VT-x
L1d cache:           32K
L1i cache:           32K
L2 cache:            256K
L3 cache:            12288K
NUMA node0 CPU(s):   0-11

Thanks.

Knowing a high-performance language will make a huge difference here. See the example below in C. I haven't tested, but it should be easy to modify for your purpose.

C++, Rust, Go, Nim and Julia can be as fast.

// Download https://raw.githubusercontent.com/lh3/minimap2/master/kseq.h
// Compile with: gcc -O2 -o myprog this-file.c -lz
#include <zlib.h>
#include <stdio.h>
#include "kseq.h"
// STEP 1: declare the type of file handler and the read() function
KSEQ_INIT(gzFile, gzread)

int main(int argc, char *argv[])
{
    gzFile fp;
    kseq_t *seq;
    if (argc == 1) {
        fprintf(stderr, "Usage: %s <in.seq>\n", argv[0]);
        return 1;
    }
    fp = gzopen(argv[1], "r"); // STEP 2: open the file handler
    seq = kseq_init(fp); // STEP 3: initialize seq
    while (kseq_read(seq) >= 0) { // STEP 4: read sequence
        int i, sum, low;
        if (seq->qual.l == 0) continue; // no quality
        if (seq->seq.l < 15 || seq->seq.l > 30) continue;
        for (i = sum = low = 0; i < seq->seq.l; ++i) {
            int qual = seq->qual.s[i] - 33; // assuming Sanger quality
            sum += qual;
            if (qual < 20) ++low;
        }
        if (sum < seq->seq.l * 30.0 || low > 2) continue;
        printf("@%s\n", seq->name.s); // output can be made faster
        puts(seq->seq.s);
        printf("+\n");
        puts(seq->qual.s);
    }
    kseq_destroy(seq); // STEP 5: destroy seq
    gzclose(fp); // STEP 6: close the file handler
    return 0;
}

To apply the program to all *.fq in a directory (again, not tested):

mkdir -p filtered
ls *.fq | xargs -i echo ./myprog {} \> filtered/{} | sh

You can also replace the final sh with parallel, but the above should be fast enough.


How to interpret Mendelian randomization results? I did Mendelian randomization using this software:
https://cran.r-project.org/web/packages/MendelianRandomization/vignettes/Vignette_MR.pdf

library(MendelianRandomization)
f=read.table("246LDout272Biobank_Retina.txt", header=T)
> head(f)
          rs exposure.beta exposure.se outcome.beta outcome.se
1  rs1029830      0.723525  0.03026430  0.066715400  0.0359278
2  rs1029832      0.723785  0.03029603  0.064105600  0.0359021
3 rs11078374     -0.411789  0.04189295 -0.000376929  0.0406439
4 rs11078382      0.882549  0.14275799 -0.197074000  0.1247720
5  rs1124961     -0.333763  0.05589377 -0.075468600  0.0576012
6  rs1135237     -0.316831  0.05552530 -0.074086200  0.0573111

MRInputObject <- mr_input(bx = f$exposure.beta,bxse = f$exposure.se,by = f$outcome.beta,byse = f$outcome.se)
EggerObject <- mr_egger(MRInputObject,robust = FALSE,penalized = FALSE,correl = FALSE,distribution = "normal",alpha = 0.05)

MR-Egger method
(variants uncorrelated, random-effect model)

Number of Variants =  246

------------------------------------------------------------------
  Method     Estimate Std Error  95% CI       p-value
  MR-Egger    0.115     0.016  0.084, 0.146   3.28e-13
 (intercept) -0.010     0.008 -0.025, 0.006   0.226
------------------------------------------------------------------
Residual Standard Error :  1.106
Heterogeneity test statistic = 298.4508 on 244 degrees of freedom,     (p-value = 0.0099)
 I^2_GX statistic: 97.2%

Let's say that my exposure is called Retina and my outcome is called Biobank. Can someone please help me interpret these results in terms of 'horizontal pleiotropy' as mentioned here: https://www.ncbi.nlm.nih.gov/pubmed/29771313

It looks like you have linkage as opposed to linkage disequilibrium, therefore horizontal pleiotropy is an explanation of the calculation.

The calculation appears to be a GLMM around the beta distribution, the residual looks low (which is good). The heterogeneity statistic is significant, you would need to assess whether this is an index of heterozygosity, but it will relate to expectation of Mendellian inhertiance based around the Hardy Weinberg equilibrium.


Restriction enzymes with multiple recognition sequences In a bioinformatics programming course I am currently teaching (using Python, and approached from a mostly CS rather than a biology background) I assigned a project which involves finding sites where restriction enzymes bind to a DNA sequence. The project is heavily influenced by Chapter 9 of the now classic text "Beginning Perl for Bioinformatics" by James Tisdall. I am using the REBASE database, in particular REBASE Format #9. In this format, each restriction enzyme has a recognition sequence which looks like e.g. C^YCGRG where ^ gives information about the exact cut-point and upper-case letters such as Y and R other than ACGT are ambiguity codes. The point of the project is to turn these into regular expressions where cleavage information is discarded and the ambiguity codes are expanded into character classes, and then use these regular expressions to search for matches.

All fairly straightforward, but one very observant student noticed that some enzymes had multiple sequences, for example M.BssHIII has a total of 5 sequences: ACGCGT,CCGCGG,RGCGCY,RCCGGY,GCGCGC. It is easy to see that these correspond to a total of 11 sequences (after the Rs and Ys are expanded), far fewer than the 144 sequences which would match VSSSSB (which is the least ambiguous single sequence capable of matching the above).

My question is, what is going on here? Am I safe in simply replacing , by | (for "or") in the regular expression, or is there something more subtle going on? Only 8 of the nearly 5000 enzymes have such multiple sequences, and neither the Tisdall book nor even the description of the format (both of which seem to assume single sequences) shed any light on the situation.

I would honestly just remove the non-homologous cutting site enzymes from the database. Firstly, these are biology students and require a simple picture to grasp the concept of coding via regex pattern matching. Secondly, in a wet-lab nobody uses these enzymes, I was not even aware the "8" existed until your post. In a teaching lab the best motto is "keep it simple".


How to release an R package with several bed files? I'm currently creating an R package, and my scripts require before any analysis is done that 1-2 bed files be loaded. Normally, I would run the scripts with the following:

library(data.table)
session_bed = read.table("/path/name/to/bedfile.bed")

and then refer to session_bed throughout the analysis.

However, in order to release this as an R package, I believe that these files must be downloaded at the time on the package installation.

(1) Normally, how are bed files distributed in R packages? Are these gzipped? Where are they placed?

According to http://r-pkgs.had.co.nz/data.html, all .bed files should be inst/extdata I guess. However, it's not clear to me how to "link" to this data within the R package.

How have other packages dealt with this? (Is there a size limit?)

(2) In this case, should the function loading the bed file be executed about loading the R library, library(package_name)? Is there more appropriate option?

You can store it in data/ subfolder. If one of your functions need this bed file, you can include importing it in your function.

Here are some examples how to get the data from data/ subfolder.

There are probably more ways to do this, this is just one way.


Dealing with indels when converting a position to chromStart and chromEnd for a bed file Starting with a chromosome # and position, I am trying to get chromStart and chromEnd values for a .bed file, but I am not sure how to calculate chromEnd when I have a variant that is an insertion or deletion.

My file is not a VCF file but it gives the chromosome number, chromosomal position (pos), reference (ref) allele, and alternate (alt) allele for each genetic variant, like you would find in a VCF file.  The tricky part is that not all of my variants are single nucleotide base changes--I also have insertions and deletions that can be multiple bases long.

For SNPs I am setting chromStart = pos - 1, and chromEnd = start + 1 = pos.  (Note, .bed format requires 0-based indexing and the end position is exclusive. See https://genome.ucsc.edu/FAQ/FAQformat.html)

I am unsure of how to handle indels though (insertions and deletions).  Currently I am adding the length of the reference allele to the starting position to get the end position (chromEnd = start + len(ref)).  I am not sure if this is right though.

I will be using this .bed file with a .bed file from Gencode and bedtools. Does anyone know the standard way the research community (or at least Gencode) treats indels when determining chromosomal position for .bed files? Thanks!

You can convert to VCF and use BEDOPS vcf2bed --insertions, vcf2bed --deletions and vcf2bed --snvs to segregate the VCF into three subsets.

For convenience, BEDOPS conversion tools deal with index shifts, including vcf2bed. The resulting BED file is 0-indexed, half-open.


Samtools Index: Chromosome Blocks not Continuous I am working with short-read whole-genome sequences from the NCBI's SRA. I have aligned and sorted all of my short-read sequences and am attempting to index each sequence into .bai format using samtools index, but am running into a couple of errors.
I unpacked the original .sra files in the following manner:
fastq-dump --defline-seq '@$sn[_$rn]/$ri' --split-files SRR6509138.sra

I then aligned each fastq paired-end read, removed duplicate reads, and converted each to bam format like so:
bwa mem xiphophorus_birchmanni_10x_12Sep2018_yDAA6.fasta SRR6509136_1.fastq SRR6509136_2.fastq | samblaster -e -r| samtools view -Sb - > blasted_SRR6509136.bam

Some of my files did not come as paired-end reads. For those files, I used the --ignoreUnmated flag on samblaster.
I then sorted each bam file like so:
samtools sort blasted_SRR6649368.bam -o sorted_SRR6649368.bam -n

I am attempting to obtain a bai file for each bam using the following command:
samtools index sorted_SRR6649368.bam sorted_SRR6649368.bam.bai

This is the error I run into for the unpaired reads:
[E::hts_idx_push] Chromosome blocks not continuous
[E::sam_index] Read 'HWI-ST387:164:D0CJWACXX:4:1101:1129:13264' with ref_name='ScdB1pO_646;HRSCAF=880', ref_length=33092979, flags=16, pos=23260108 cannot be indexed
samtools index: failed to create index for "sorted_SRR791885.bam": No such file or directory

This is the error I run into for the paired reads:
[E::hts_idx_push] Unsorted positions on sequence #79: 11159020 followed by 11158717
samtools index: failed to create index for "sorted_SRR6649368.bam"

Can anyone help me figure out why this is happening?

samtools sort blasted_SRR6649368.bam -o sorted_SRR6649368.bam -n

These error messages indicate that the reads are not sorted by coordinate — in particular, that the reads mapped to ScdB1pO_646;HRSCAF=880 are not all together, and that the reads at positions 11159020 < 11158717 are not sorted by position on their chromosome.
This is because samtools sort -n has been used to sort the reads by name instead. Remove -n to sort by position, which is what is needed to prepare a BAM file for indexing with samtools index.


OpenAPI-Generator Python Client Instructions Examples Trying to learn a little bit of Python and Biopython I ended up figuring out what a Rest API is when the RCSB PDB repository released its new Rest API.
I tried to retrieve info from the PDB repository using python3 requests,
then I managed to use OpenAPI-generator to generate what they call a 'Python3 Client'.
Now I am at loss I have this ApiClient class:
print dir():

 ['NATIVE_TYPES_MAPPING', 'PRIMITIVE_TYPES', '_ApiClient__call_api',
  '_ApiClient__deserialize', '_ApiClient__deserialize_date',
  '_ApiClient__deserialize_datetime', '_ApiClient__deserialize_file',
  '_ApiClient__deserialize_model', '_ApiClient__deserialize_object',
  '_ApiClient__deserialize_primitive', '__class__', '__delattr__',
  '__dict__', '__dir__', '__doc__', '__enter__', '__eq__', '__exit__',
  '__format__', '__ge__', '__getattribute__', '__gt__', '__hash__',
  '__init__', '__init_subclass__', '__le__', '__lt__', '__module__',
  '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__',
  '__setattr__', '__sizeof__', '__str__', '__subclasshook__',
  '__weakref__', '_pool', 'call_api', 'client_side_validation',
  'close', 'configuration', 'cookie', 'default_headers',
  'deserialize', 'files_parameters', 'parameters_to_tuples',
  'pool', 'pool_threads', 'request', 'rest_client',
  'sanitize_for_serialization', 'select_header_accept',
  'select_header_content_type', 'set_default_header',
  'update_params_for_auth', 'user_agent']

and don't know how to use it. Googling for examples or searching for info on the OpenAPI-generator website doesn't seem to work for me, I find plenty of tutorials and examples on how to generate the client but none about how to use it.
The generated client class should help to better interact with the server as I understood (I could be wrong here). I am asking for helps and hints here because the OpenAPI webpage under Community points to stackoverflow under [openapi-generator] tag. Any tutorial, easy documentation that could help me understand the output of the generator?

OK fine, found useful information on READme.md generated by OpenAPI-generator itself
and tried pdoc to get html document of my OpenAPI-generator generated client Module
Moreover on stackoverflow How to generate API Documentation using openapi-generator? you can find how OpenApi Generator can do generation of html documentation from your yaml or json definition of your API.
For openAPI json api specification I found RapiPDF too : Web Component to generate PDF from OpenAPI Spec (json yml file).


Calculating polarity and hydrophobicity of a ligand/small molecule I have lot of SMILES IDs of the ligands. What I am interested in is the polarity and hydrophobicity score of each of the ligands.

In other words, does ligand have more polar groups and can it make more polar interactions or is it more hydrophobic?

Does anybody know the way/software to look at it?

Best

You can calculate their polar surface area or logo using rdkit or openbabel. Other molecular properties may be available as well.


Filtering VEP annotation file I am filtering a VEP annotated vcf, trying to maintain just those variants classified as deleterious by SIFTand as damaging (probably or possibly included) by PolyPhen.

I am using:

filter_vep -i "$input" -o "$output" -filter "SIFT = deleterious and PolyPhen match damaging"

and I also tried:

filter_vep -i "$input" -o "$output" -filter "SIFT match deleterious and PolyPhen match damaging"

The output contains just variants with some value in these two fields, but some of the variants are classified as tolerated by SIFT and as benign by PolyPhen, so do you know why the filter is not working?

As an example: this two variants shouldn't exist in my filtered VCF, because the first one is classified as benign by PolyPhen and the second is classified as benign by both fields:

    CHROM  POS  REF  ALT  FILTER  INFO
    1       1934441 .       G       T       .       PASS    deleterious(0.01)|benign(0.138)
    1       207085100       .       G       A       .       PASS   tolerated(0.16)|benign(0.037)

As an example of a variant that has pass properly those filters:

    CHROM  POS  REF  ALT  FILTER  INFO
    2       3938449 .       T       A       .       PASS    deleterious(0.01)|probably damaging(0.67)

Those are 3 variants among those that I have in my filtered VCF, but I should be getting just those variants which are classified as deleterious and damaging by both SIFT and PolyPhen, as the third variant I've shown.

If you are working under a Linux environment, then you may easily use a piped grep search. For example:

grep deleterious file.vcf |grep damaging >filtered_file.vcf


Why do I get so many insertions from Minimap2 on my Nanopore WGS? I'm a starting my analysis on nanopore whole genome sequencing. I start my analysis from this popular Github.

The sample I downloaded was WGS for NA12878, so I would assume it's alignment to GRCh38 shouldn't be that bad.

But ... I'm getting lot's of insertions...

I'm not sure it was caused by the sequencing platform, the MiniMap aligner or my bugs.

I supply all my commands, starting from sample download.

wget https://s3.amazonaws.com/nanopore-human-wgs/rel5-guppy-0.3.0-chunk10k.fastq.gz
gunzip rel5-guppy-0.3.0-chunk10k.fastq.gz
wget http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa
minimap2 -d GRCh38_full_analysis_set_plus_decoy_hla.fa.index GRCh38_full_analysis_set_plus_decoy_hla.fa
minimap2 -ax map-ont -t 20 GRCh38_full_analysis_set_plus_decoy_hla.fa.index rel5-guppy-0.3.0-chunk10k.fastq | samtools view -bS - | samtools sort -@ 5 > rel5-guppy-0.3.0-chunk10k.bam
samtools index rel5-guppy-0.3.0-chunk10k.bam

Insertions and deletions are the dominant error mode of long read sequencing, including nanopore sequencing. What you see is not unexpected. Things may have improved by now if you would download the raw fast5 data and repeat the basecalling. There is no need to gunzip the fastq.gz prior to alignment. Your commands for alignment look alright, except that (if you are using a recent version of samtools) you don't really need samtools view and can pipe directly to samtools sort, using -o to specify the output format:

minimap2 -ax map-ont -t 20 <index> <fastq> | samtools sort -@5 -o alignment.bam

In a recent version of igv you can right-click the alignment and chose a "Quick consensus mode" (or something like that), which will make your alignments at least a bit prettier. Note how the consensus sequence (top of the graph) is quite clean - while individual reads may be noisy this is mostly randomly distributed and as such ironed out in consensus.


Do any computational phylogenetic methods enable the specification of ancestral states? Various phylogenetic algorithms estimate ancestral states of a phylogenetic dataset. Is there a way in either maximum parsimony, distance-based methods, or Bayesian inference to indicate what the ancestral states of characters were?

For maximum parsimony, Fitch's method seems simple and intuitive, and can be done in the R phangorn package: Fitch, Walter M. (1971). "Toward Defining the Course of Evolution: Minimum Change for a Specific Tree Topology". Systematic Zoology. 20 (4): 406. doi:10.2307/2412116. ISSN 0039-7989. JSTOR 2412116

There are various methodological criticisms of parsimony that you can find in the literature, but its appeal is you can easily understand why the algorithm came to the conclusions it did.

RAxML can do maximum likelihood reconstructions based on a tree and tip values, and I think Mr. Bayes can do a Bayesian version.


snakemake with --conda-prefix, can I use a pre-built conda environment? I have read this question (Running Snakemake in one single conda env) but still have some doubt.

I am using the '--conda-prefix /some/dir' option and I have a rule in my snakemake file with the directive:

rule xxx:
  ...
  conda:
    "envs/some.yaml"
  ...

The first call of that rule will create the conda environment according to the yaml file in /some/dir/envs with some random name (in my case 'bf55a79e'). If I run the workflow in a different directory, the same conda environment (/some/dir/envs/bf55a79e) gets activated directly instead of being built again. So this is what I expected.

My question is:

is there some way of skipping the building and specifying that rule xxx should use an already existing conda environment somewhere?
does anyone know where the mapping between rule and directory gets stored?

If the environment already exists it will simply be used. You cannot, however, tell it to use an environment with a normal name (e.g., "my_env").
I assume you mean, "where the heck does the random looking conda environment name come from?!?", which is a very good question. The answer to that is that it's a hash of the --conda-prefix setting and the contents of the environment yaml file you specify. The actual code for doing this (I'm "stealing" it from snakePipes) is:

import hashlib
md5hash = hashlib.md5()
md5hash.update("what you gave to --conda-prefix".encode())
f = open("your environment yaml file", 'rb')
md5hash.update(f.read())
f.close()
h = md5hash.hexdigest()

h is then the directory name. There are actually two lengths of those, the shorter 8 character version that snakeMake uses by default or the full length version that the code above will produce. snakeMake will actually use either of these. If you're wondering, no I don't think any of this is documented anywhere, you have to dig through the snakeMake code to find it.


Multi-line efetch digging I'm using efetch to dig out a middle "envelope" protein within a viral genome specifically within fasta format. The code below works but may not be robust for an entire database. Thus I want to replace the two "perl" lines with formal code within the Entrez package. I tried using efilter -query but couldn't get it to work. Genome positions are one approach but are not absolute from one virus to another. Rewriting this may be the best option.

ncbi2=NC_002021
efetch -db protein -format gpc -id "$ncbi2" |\
xtract -pattern INSDSeq -element INSDQualifier_value |\
perl -p -e 's/^(.+?)\s+taxon.*.*\s+\b([A-Z]{300}[A-Z]+).*/>$1\n$2/'|\
perl -p -e "s/[\t \;\/]+/_/g"

Output
>InfluenzaAvirus_A_PuertoRico_8_1934
MDVNPTLLFLKVPAQNAISTTFPYTGD.....

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

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

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

For more information, checkout EDirect documentation and some  example queries


Plot percentage of genome covered Given a aligned bam file (wgs/bwa-mem) what steps do i need to perform to generate a plot as seen below (from this paper):

We are trying to see how much of genome was covered using pacbio and illumina for cow.

You mean remove this left plot?

The general steps are:

Install deepTools
Run plotCoverage on your BAM files
Remove the left-most plot in photoshop/gimp/imagemagick


Generating multiple NCBI URLs with eutils I would like to create a URL that takes input from different arrays to generate a search. The data I would like to retrieve is from NCBI. Below you can find the code I wrote. But I can't get it to work; I am stuck. Thank you in advance.

I should have something like this as result but i got nothing. https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=gene&retmode=xml&retmax=10&sort=relevance&term=A20%20AND%20Homo%20sapiens%20%5borgn%5d%20AND%20alive%5bprop%5d

$eutils = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils";
$db = "gene";
$gene_symbols = array("A20","HBB","ABO");
$species = array("homo sapiens", "Escherichia coli","Saccharomyces cerevisiae");
$links = array();
foreach ($gene_symbols as $key => $GS){
    $SP = $species[$key];
    $query = "$GS AND $SP [orgn] AND alive";
    $esearch = "$eutils/esearch.fcgi?db=$db&retmode=xml&retmax=10&sort=relevance&term=";
    #echo $esearch;
    $links[] = get($esearch.$query);
}
foreach($links as $lk){
    xml = simplexml_load_file($lk) or die("feed not loading");
}
var_dump(xml);

As output I got nothing.

But I would like to have this: https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=gene&retmode=xml&retmax=10&sort=relevance&term=A20%20AND%20Homo%20sapiens%20%5borgn%5d%20AND%20alive%5bprop%5d
And this for every gene. For some reason it put them all after each other and the foreach loop doesn't seem to work. I got no errors.

I run this code in php using visual studio code. The result should be visualised in the webbrowser (http://localhost/test.php).

I don't know anything about php, but when I copy your code into here it isn't happy about the variable name $query.

This code appears to produce valid urls, but I can't make the get part work:

$eutils = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils"; $db = "gene";$gene_symbols = array("A20","HBB","ABO");
$species = array("homo sapiens", "Escherichia coli","Saccharomyces cerevisiae"); $links = array();

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

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

   echo $esearch.$SearchTerm;
}

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


Get regions' information from DNA sequence data (bsgenome.hsapiens.ucsc.hg19) I have a problem in R. I have the following dataSet (the first three rows shown) (the 5th number is the methylation level in its region):

chr1    15864   15866   FALSE   894 +
chr1    534241  534243  FALSE   921 -
chr1    710096  710098  FALSE   729 +

and then I'm trying to get the sequences from every row in the dataSet (427.000 rows in total):

library(BSgenome.Hsapiens.UCSC.hg19)
genome <- BSgenome.Hsapiens.UCSC.hg19
chr<-as.matrix(as.character(mydataSet[,1]))
start<-as.matrix(as.integer(as.character(mydataSet[,2]))-50)
end<-as.matrix(as.integer(as.character(mydataSet[,3]))+50)
Sequences_for_50_plus_bases<-getSeq(genome,chr,start=start,end=end)

I'd like to find a new fuction where I can have this dataSet, and find for its rows (which represent sequences in the human genome) if they refer to a promoter, or an enhancer or any other DNA regions' information, because I'd like to compare the differences between the regions' methylation levels.

Is there a way to do this?

You can get specific sequences efficiently using the following packages:

library(GenomicFeatures)
library(Biostrings)
library(tidyverse)
library(data.table)
library(dtplyr)
library(BSgenome.Hsapiens.UCSC.hg19)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)

Methylation_data <- tribble(
  ~seqnames,~start,~end,~Logical,~"Methylation Level",~strand,
  "chr1",15864,15866,FALSE,894,"+",
  "chr1",534241,534243,FALSE,921,"-",
  "chr1",710096,710098,FALSE,729,"+") # Arrange your dataset

Methylation_GRanges <- GRanges(  # Construct a GRanges objects
  seqnames = Methylation_data %>% pull(seqnames), # First pull in the chromosome locations
  IRanges(start = (Methylation_data %>% pull(start))-50, # then the start position minus 50 nt
          end = (Methylation_data %>% pull(end))+50), # the end position plus 50 nt
  strand = Methylation_data %>% pull(strand) # the strand your sequence is located on
)
> getSeq(Hsapiens,Methylation_GRanges)
A DNAStringSet instance of length 3
width seq
[1]   103 GCTGCTGCTTCTCCAGCTTTCGCTCCTTCATGCTGCGCAGC...CTTGGCGGATGGACTCTAGCAGAGTGGCCAGCCACCGGAG
[2]   103 ACAGTTCCAATGTAATCAGAGAGAACATCACACACACACCA...ACCTGAGCAGCACTCTGCAAAGCTGTCAAGGCGGTGAAGC
[3]   103 GATAAAATAAAGCTTAGATTGGAAAAAATATTTAAGATTCT...GGAAGCTGAGTAATTGTATGTTCAAATACTTGCAAAACAT

To write the sequences to a file, use writeXStringSet()

For the second part of what you need, GenomicFeatures has a number of functions that obtain the annotated sequences of the whole genome by its type:

transcriptsBy(x, by=c("gene", "exon", "cds"), ...)
exonsBy(x, by=c("tx", "gene"), ...)
cdsBy(x, by=c("tx", "gene"), ...)
intronsByTranscript(x, ...)
fiveUTRsByTranscript(x, ...)
threeUTRsByTranscript(x, ...)

And this is how you use one of them to find whether your sequences overlap a particular type of DNA:

Find_overlap <- function(Ref_Data,
                         Overlap_Message){

  Ref_Data <- as.data.table(Ref_Data) # Convert to data.table for efficiency

  Methylation_data <- tribble(
    ~seqnames,~start,~end,~Logical,~"Methylation Level",~strand,
    "chr1",15864,15866,FALSE,894,"+",
    "chr1",534241,534243,FALSE,921,"-",
    "chr1",710096,710098,FALSE,729,"+") %>%
    mutate(Overlap = "-",
           start = start-50,
           end = end+50) %>%
    as.data.table()

  Methylation_data[Ref_Data,
                   Overlap := Overlap_Message,
                   on = .(start >= start,
                          end <=  end)] # wherever the range overlaps, assign an identifier to it
}

Methylated_in_5UTR <- fiveUTRsByTranscript(TxDb.Hsapiens.UCSC.hg19.knownGene) %>%
  Find_overlap(., "5' UTR") %>%
  print

seqnames  start    end Logical Methylation Level strand Overlap
1:     chr1  15814  15916   FALSE               894      +       -
2:     chr1 534191 534293   FALSE               921      -       -
3:     chr1 710046 710148   FALSE               729      +  5' UTR
```


residual Squared Coefficient of Variation (rCV²) vs Distance to Median (DM) I am working with single cell RNA-seq data. I obtained the squared coefficient of variation (CV²) as a measure of gene expression variability:

I want a metric to express gene expression noise that eliminates/takes into account the observed relationship between CV² of a gene and its mean expression so I can compare this metric across genes/conditions:

So far I have found two metrics that try to tackle this problem:

Distance to median (DM):

  A median-based trend is fitted to the log-transformed CV² against the log-transformed mean. The DM is defined as the residual from the trend for each gene. This statistic is a measure of the relative variability of each gene, after accounting for the empirical mean-variance relationship

residual CV² (rCV²):

  Fit between the mean log2 expression and CV² using a gamma generalised linear model as parametrisation then computing the absolute deviation of the CV² values from the fitted CV² values.

Do you have any clue about which metric is better or if there exists any other similar metrics?

There has been a lot of work done on this problem already, particularly these two papers:

https://www.biorxiv.org/content/10.1101/576827v2

https://www.biorxiv.org/content/10.1101/574574v1


Is there an "official" standard usage of Red and Green on differential gene expression heatmaps? I have recently found myself making multiple heatmaps for visualization of differential gene expression results between replicates in two experimental conditions.
The default settings in the plotting software I use color:

Upregulated genes --> red
Downregulated genes --> green
Unchanged genes --> black

My intuition screams at me that this color scheme is reversed, and that upregulated genes should be green and that downregulated genes should be red. I have asked several labmates about their intuition regarding this color scheme, and have found about a 50/50 split in those that intuitively agree with each different scheme.
In my searching on this topic, I have found that several programs plot heatmaps with the same default colors as the program I'm using (red=upregulated, green=downregulated). Additionally, I have been told in conversation that this is essentially an accepted standard.
If this is indeed the convention, I will begin trying to rewire my intuition, but I want to confirm that it is before I do so. However I have not yet been able to find any "official" guidelines from agencies, publishers, authorities in the field, etc., specifying the use of red/green in heatmaps.
I do know that there are many suggestions for alternative heatmap color schemes, but I'm interested in finding an official answer for the red/green scheme specifically.
Edit: Several comments have stated that there is no such standard as described above. Therefore I am modifying the following question somewhat to:
Are there "official" guidelines from agencies, publishers, authorities in the field, etc., specifying the use of red/green in heatmaps?
Or conversely are there any “official” guidelines stating that there is no such specification and/or that any such selection of colors is acceptable?

As others already pointed out, there is no official heatmap color scheme. I am not sure there is an official color scheme for any other types of plots either.
However, heatmaps seem to be far more standardized. Traditionally, heatmaps have been red-black-green as you originally stated. Most heatmaps from the microarray era and early RNA-seq are like that. You can have a heatmap for any type of data, but the majority is for gene expression. A lot of tools and tutorials were written back then, so they have that color scheme. From about 2008 to 2012, most people switched to the red-white-blue (or red-yellow-blue) color scheme. Switching green to blue avoids the color-blindness issue. Switching black to white is better for the ink if you print your papers and generally looks better on white background.
I originally found red as high to be confusing because I was thinking in financial terms (if the market is red, it's down). Now I think of it as a water faucet (red is hot, blue is cold).
Edit: In the last year or so, the viridis palette started to become widely used. It's far from the standard, but could become more popular in the future.


Is it possible to perform MinION sequencing offline? I vaguely remember, that the original plan of Oxford Nanopore was to provide cheap sequencers (MinION), but charge for base-calling. For that reason the base-calling was performed in the cloud, and the plan was to make it commercial once the technology is established.

Of course, this limits the potential uses of MinION in the field, since huge areas do not have decent internet connection. Also, not all the data can be legally transferred to a third-party company in the clinical studies.

For example for the Ebola paper, they had to create a special version of their software:

  An offline-capable version of MinKNOW, with internet ‘ping’ disabled
  and online updates disabled was made available to us by Oxford
  Nanopore Technologies specifically for the project

There are couple of third-party base-callers available today. I am aware of Nanocall and DeepNano, but since they are not official, it can be hard for them to keep-up with the latest versions of sequencers and cells.

Is it possible as of today to sequence offline without a special arrangement (like the Ebola one).
If not, what's the policy of Oxford Nanopore toward third-party base-callers? Are they going to help them, or to sue them eventually?

Short answer: yes, but you need to get permission (and modified software) from ONT before doing that.

... but that doesn't tell the whole story. This question has the potential to be very confusing, and that's through no fault of the questioner. The issue is that for the MinION, sequencing (or more specifically, generating the raw data in the form of an electrical signal trace) is distinct and separable from base calling. Many other sequencers also have distinct raw data and base-calling phases, but they're not democratised to the degree they are on the MinION.

The "sequencing" part of MinION sequencing is carried out by ONT software, namely MinKNOW. As explained to me during PoreCampAU 2017, when the MinION is initially plugged into a computer it is missing the firmware necessary to carry out the sequencing. The most recent version of this firmware is usually downloaded at the start of a sequencing run by sending a request to ONT servers. In the usual case, you can't do sequencing without being able to access those servers, and you can't do sequencing without ONT knowing about it. However, ONT acknowledge that there are people out there who won't have Internet access when sequencing (e.g. sequencing Ebola in Africa, or metagenomic sequencing in the middle of the ocean), and an email to <support@nanoporetech.com> with reasons is likely to result in a quick software fix to the local sequencing problem.

Once the raw signals are acquired, the "base-calling" part of MinION sequencing can be done anywhere. The ONT-maintained basecaller is Albacore, and this will get the first model updates whenever the sequencing technology is changed (which happens a lot). Albacore is a local basecaller which can be obtained from ONT by browsing through their community pages (available to anyone who has a MinION); ONT switched to only allowing people to do basecalling locally in about April 2017, after establishing that using AWS servers was just too expensive. Albacore is open source and free-as-in-beer, but has a restrictive licensing agreement which limits the distribution (and modification) of the program. However, Albacore is not the only available basecaller. ONT provide a FOSS basecaller called nanonet. It's a little bit behind Albacore on technology, but ONT have said that all useful Albacore changes will eventually propagate through to nanonet. There is another non-ONT basecaller that I'm aware of which uses a neural network for basecalling: deepnano. Other basecallers exist, each varying distances away technology-wise, and I expect that more will appear in the future as the technology stabilises and more change-resistant computer scientists get in on the act.

Edit: ONT has just pulled back the curtain on their basecalling software; all the repositories that I've looked at so far (except for the Cliveome) have been released under the Mozilla Public License (free and open source, with some conditions and limitations). Included in that software repository is Scrappie, which is their testing / bleeding-edge version of Albacore.


Relationship between sequencing lane and ngs dataset I am fairly new to NGS data analysis and I am struggling to understand the exact relationship between a sequencing lane and an NGS dataset. I should add I don't work in the lab, I only do bioinformatics.
I understand the basics of how to make NGS libraries (very briefly: shearing, ligation of adapters, optional PCR, ligation to solid surface, following by sequencing-by-synthesis or other method). What I am struggling with is the visualisation of what happens once you have made your libraries. My questions are:

If you have used adapters containing barcode sequences (that identify each patient), you would presumably end up with one single tube with your libraries in it right?

If so, you then "put" the contents of that tube along with DNA polymerase and nucleotides into one single lane of a flow cell? Your provider would then give you raw FASTQ files that have been identified using this barcode? So you would effectively received one (or 2, depending if it's single/paired-end sequencing) FASTQ file per patient right?

If you do not use adapters with barcodes, does that mean you would have to use a sequencing lane per patient? If so, that would vastly increase the price you pay your provider, would it not?
Apologies if these are basic questions.

If you have used adapters containing barcode sequences (that identify each patient), you would presumably end up with one single tube with your libraries in it right?

Correct, you have to prep each sample into a library.

  If so, you then "put" the contents of that tube along with DNA polymerase and nucleotides into one single lane of a flow cell?

You could, or you could pool it with other samples (assuming they have compatible indices). As @swbarnes2 mentioned, you don't see a lot of samples run on a single lane these days. Sequencing centers like to split things across lanes in case on lane fails for some reason or there are bubbles in the flow cell or other technical issues.

  Your provider would then give you raw FASTQ files that have been identified using this barcode?

Typically, yes, and if your sample was run on multiple lanes they'll often be merged for you automatically. This isn't always the case. I've seen a number of people receive samples either still split by lane (annoying, but you can just cat things together) or, much worse, not demultiplexed at all. In the latter case, people often receive 3 (or 4) fastq files per lane, one for each read and another for each index read. It's rare that sequencing facilities do this, but it happens.

  If you do not use adapters with barcodes, does that mean you would have to use a sequencing lane per patient? If so, that would vastly increase the price you pay your provider, would it not?

You could sample a sequence per lane. Whether that makes economic sense is up to you to decide.


What is mate rescue in bwa mem? BWA mem has the -S and -P tags for skip mate rescue and skip pairing; mate rescue performed unless -S also in use.

What do these do? I presume -P aligns read pairs independently of each other. Is that correct?

And what does -S do?

Using -SP is equivalent to running bwa mem on each of the two mates as if they are single-end reads, but it formats the output as a proper paired-end output, i.e. with all pair-related flags added properly. Without -SP, by default bwa mem forces an alignment of a poorly aligned read if its mate is aligned somewhere. -SP turns off the forced alignment. We use -SP for mapping Hi-C reads.


What pitfalls exist with running FastQC on a bam file? I was analysing a bad sequencing run of some RNA data using FastQC, I supplied it RNA-STAR-aligned bam files (and the human reference).  The output of MultiQC had many more "unique reads" in the sequence counts table than there were reads in the FastQ files that were fed to RNA-STAR.
At first, I realized that all the FastQC stats were about "mappings", not reads (as they were labeled), but "unique mappings" doesn't seem to make sense either, because even if I consider the possibility that segments of reads were mapping in different places, the numbers don't make sense.  Here is the worst case:

Max read length after trimming: 47 bases
Number of reads in the original FastQ file that went into STAR: 59M

Bam Stat Table Numbers:

Total records: 391 million
Unique mappings: 3.1 million
Non-primary Hits: 335 million

FastQC Sequence Counts Bar graph:

"Unique reads": 185 million
"Duplicate reads": 206 million

While I found documentation stating that you can run FastQC on bam files, there were no caveats I could find about it needing to be unmapped reads.
Can anyone explain precisely why there can be 185 million "unique reads" when bam stats says there were 3 million unique mappings and the number of reads is 59M?

I don't see why you are trying to compare two totally different measurements.
When FASTQC says a read is unique, it means that no other read shares its sequence.
When STAR says a read is uniquely mapped, that means it maps to one and only place.  That does not mean that there aren't a hundred more identical reads that map nicely to the exact same place.
You could have a very high level of PCR duplication of reads that nicely align to one and only one place. Or you could have a whole bunch of reads that differ slightly from each other who all map to 5 locations each.
FASTQC was designed for DNA sequencing, not RNASeq.  Some of its warnings make no sense in the context of RNASeq.  If I was doing RNASeq of a filtered set of pancreatic alpha cells,  I would expect a hell of a lot of reads for glucagon, and a lot of them would be natural duplicates, and that would be biologically normal, and not a problem.
You could take the same batch of reads, and align them to one genome, and align them to a mix of two closely related genomes.  The second will have far more non-uniquely mapped reads, but the FASTQC stats will not have changed.
That said, if you doing human or mouse RNASeq, 1% unique alignment means something is wrong.  Find out where all those reads are aligning, you might have rRNA contamination.  If you are doing an organism with a high level of duplication in its genome that might be normal.


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

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


Is there a standard k-mer count file format? I am doing a research project involving calculating k-mer frequencies and I am wondering if there is any standard file format for storing k-mer counts.

Not as far as I am aware. The Ray assembler used to (and possibly still does) store the kmers as FASTA files where the header was the count of the sequence, which I thought was a pretty neat bastardisation of the FASTA file format. It looks like this format is also used by Jellyfish when reporting kmer frequencies by the dump command (but its default output format is a custom binary format):

  The dump subcommand outputs a list of all the k-mers in the file associated with their count. By default, the output is in FASTA format, where the header line contains the count of the k-mer and the sequence part is the sequence of the k-mer. This format has the advantage that the output contains the sequence of k-mers and can be directly fed into another program expecting the very common FASTA format. A more convenient column format (for human beings) is selected with the -c  switch.

Jellyfish changed their internal format between v1 and v2 (both not FASTA), because they changed to doing counts based on bloom filters. Jellyfish2 has an optional two-pass method that sets up a bloom filter intermediate file to record kmers, and multiple different final reporting formats.

Khmer also uses bloom filters, but in a slightly different way. It also has been extended to be useful for partitioning and comparing datasets.


Importance of Proper Pairs vs Aligned Reads for RNASeq data I have stranded, paired end RNASeq reads that I have aligned using STAR. I plan do conduct a differential expression analysis with DESeq2. After running quality control checks, a good portion of my reads are aligned but I'm a little concerned about the proportion of properly paired reads

Here is the inner distance plot for the samples

For human RNA sequencing, does anyone have a rule of thumb that they'd consider is an appropriate proportion of properly paired reads? Are there any good papers that talk about this issue? As MultiQC marks the column as yellow, I'm concerned the proper pair percentage may not be satisfactory.

How important are properly paired reads for differential expression analyses?

If you check the RSeQC docs you can get an explanation of how the inner distance stats work. By the look of your inner-distance MultiQC plot I guess that your sequencing is 2x100bp reads. This plot shows that your insert sizes on the pairs is pretty short, meaning that a lot of the paired reads are overlapping significantly. This will be affecting the alignments a little, resulting in the low proportion of properly paired reads.

The good news is that for differential expression, this probably doesn't really matter too much. What is important is getting fragment counts, so as long as at least one of the pair mates is mapping properly, then you will will be getting sensible counts - even if read 2 is not mapping properly. You should be able to have confidence in your results from multiple replicates if in doubt - you have more than just -rep1 right? ;)

Note that the yellow background in the MultiQC column is not passing judgement with categories like pass/warn/fail in FastQC - it just happens that the colour scheme ended up there for your samples. It's a continuous gradient from red to yellow to green for those values (0-100%).


Plotting Ramachandran Plot from more than one PDB file I want to plot a Ramachandran Plot for around 5000 PDB files in python that is saved in a folder in my system. I am searching for a way to make it without accessing each file by hard-coding each file name. Please tell if there is a way around this.

You could probably use PyRAMA for that. The only problem is that by PyRAMA plots to display, while you probably would like to save to a file.

One way to get around this is to:

Patch PyRAMA to allow ramachandran function to accept output image filename.

def ramachandran(file_name_list, output_image=None):
    ...
    plt.tight_layout()
    if output_image is not None:
        plt.savefig(output_file, dpi=300)
    else:
        plt.show()

If this works as expected I would really encourage you to submit a pull-request (or at least open an issue) on github.
Now you can write your main script:

from glob import iglob

from rama_lib.ramachandran_calc import ramachandran

for pdb_filename in iglob('/my/dir/with/pdb/', '*.pdb'):
    img_filename = pdb_filename.rsplit('.', 1)[0] + '.png'
    ramachandran((pdb_filename,), img_filename)

I did not test the code, but I hope you get an idea.

There are also other ways to make a ramachandran plot in python. Please check: 1, 2, 3.


Given a VCF of a human genome, how do I assess the quality against known SNVs? I'm looking for tools to check the quality of a VCF I have of a human genome. I would like to check the VCF against publicly known variants across other human genomes, e.g. how many SNPs are already in public databases, whether insertions/deletions are at known positions, insertion/deletion length distribution, other SNVs/SVs, etc.? I suspect that there are resources from previous projects to check for known SNPs and InDels by human subpopulations.

What resources exist for this, and how do I do it?

To achieve (at least some of) your goals, I would recommend the Variant Effect Predictor (VEP). It is a flexible tool that provides several types of annotations on an input .vcf file.  I agree that ExAC is the de facto gold standard catalog for human genetic variation in coding regions.  To see the frequency distribution of variants by global subpopulation make sure "ExAC allele frequencies" is checked in addition to the 1000 genomes.

Output in the web-browser:

If you download the annotated .vcf, frequencies will be in the INFO field:

##INFO=<ID=CSQ,Number=.,Type=String,Description="Consequence annotations from Ensembl VEP. Format: Allele|Consequence|IMPACT|SYMBOL|Gene|Feature_type|Feature|BIOTYPE|EXON|INTRON|HGVSc|HGVSp|cDNA_position|CDS_position|Protein_position|Amino_acids|Codons|Existing_variation|DISTANCE|STRAND|FLAGS|SYMBOL_SOURCE|HGNC_ID|TSL|SIFT|PolyPhen|AF|AFR_AF|AMR_AF|EAS_AF|EUR_AF|SAS_AF|AA_AF|EA_AF|ExAC_AF|ExAC_Adj_AF|ExAC_AFR_AF|ExAC_AMR_AF|ExAC_EAS_AF|ExAC_FIN_AF|ExAC_NFE_AF|ExAC_OTH_AF|ExAC_SAS_AF|CLIN_SIG|SOMATIC|PHENO|MOTIF_NAME|MOTIF_POS|HIGH_INF_POS|MOTIF_SCORE_CHANGE

The previously mentioned Annovar can also annotate with ExAC allele frequencies.  Finally, should mention the newest whole-genome resource, gnomAD.


How to output only the overlapping region of two files? I have two files of genomic coordinates. I am trying to output only the overlapping region of fileA if it intersects with fileB. Not the original fileA coordinates.

fileA

chr1 25 50
chr1 75 200

fileB

chr1 10 60
chr1 80 90

desired output

chr1 25 50
chr1 80 90

bedtools is my go-to program for operations on genomic intervals. In particular, the bedtools intersect operation is what you're looking for here.

$ cat fileA
chr1    25      50
chr1    75      200
$ cat fileB
chr1    10      60
chr1    80      90
$ bedtools intersect -a fileA -b fileB
chr1    25      50
chr1    80      90

Make sure that the values are separated by tabs rather than spaces, or bedtools will have trouble reading the files.


less than 80% of your list has mapped to your chosen identifier type I am trying to do GO analysis. However, when I run DAVID, I am getting this error:

You are either not sure which identifier type your list contains, or less than 80% of your list has mapped to your chosen identifier type. Please use the Gene Conversion Tool to determine the identifier type.
In addition, when I am trying to add a column with gene names next to esemble ids in R in the file with differential expression genes, many genes say NA under gene name column.

I have used STAR and Salmon for the alignment,and DESEQ2 for differential expression analysis.
I have used premade indexes for both tools.
Thank you in advance!
Some of the ENSEMBL codes I used:

ENSG00000249200 ENSG00000239920 ENSG00000242295 ENSG00000235763 ENSG00000226933 ENSG00000283632 ENSG00000235838 ENSG00000232429 ENSG00000281076 ENSG00000271042 ENSG00000229119 ENSG00000253384 ENSG00000285645 ENSG00000240627 ENSG00000213590 ENSG00000270739 ENSG00000240201 ENSG00000203531 ENSG00000271626 ENSG00000250316 ENSG00000218713 ENSG00000231548 ENSG00000231049 ENSG00000286185 ENSG00000214535 ENSG00000212961 ENSG00000219747 ENSG00000259100 ENSG00000231767 ENSG00000243175 ENSG00000236804 ENSG00000240634 ENSG00000270773 ENSG00000230267

DAVID last received an update in May 2016 (https://david.ncifcrf.gov/content.jsp?file=update.html), so any additions to the Gene Ontology database since then won't be present. I checked the examples you posted in Biomart and found that only 2 have GO annotations, so you're likely not missing an immense number of genes. Having said that, you'd be best served to use a more frequently updated tool than DAVID. Around 2015 there was an outcry because it hadn't been updated in 6 years and the database was so out of date that people were getting questionable results. Have a look at one of the many Bioconductor packages for convenient free solutions that are better maintained.


Find overlap between VCF files I have two VCF files and I want to compare the missing rate in each of them, but I want to only look at sites that are present in both files.  How could I go about getting a list of positions that are included in both VCF files and then filtering those VCF files to only include those positions?

Something like vcf-isec? Quoting:

  Creates intersections and complements of two or more VCF files. Given multiple VCF files, it can output the list of positions which are shared by at least N files, at most N files, exactly N files, etc. The first example below outputs positions shared by at least two files and the second outputs positions present in the files A but absent from files B and C.


Looking for a tool to find 16S rRNA in hundreds of genomes I have hundreds of bacterial genomes and I would like to know which is the best informatic tool in which I can run all genomes on one go and get a 16S rRNA gene sequence of each genome as output.

Could anyone give me some advice?

You can use Barrnap

barrnap -k bac --threads "$(nproc)" \
        -o output_rrna.fna < input_file.fna \
        > output_rrna.gff 2> /dev/null

Put the file in one line of sequence per fasta id and get sequence with 16S annotation.

perl -pe 'if(/\>/){s/\n/\t/}; s/\n//; s/\>/\n\>/' output_rrna.fa | \
    perl -pe 's/\t/\n/' | tail -n+2 | cut -d\: -f1 | \
    grep -A1 16S > output_16S.fasta

Check the output_rrna.gff and make a BLAST of the sequence to be sure about the predicted 16S.


Why does cutadapt remove low quality bases from the ends of reads only? I use cutadapt to remove low quality bases from my Illumina reads. The algorithm only removes low quality bases from the end until it reaches a good quality base. If there is a bad quality base beyond that, it is not trimmed. Why? Why doesn't the algorithm remove all low quality bases? It could replace low quality bases in the middle of a read by an N, for example. Or is it that the low quality base is still most probably correct and hence one does not want to lose the base information for mapping?

If you check the read QC statistics of an Illumina run in e.g. fastQC, you will see that at the end of the read the quality decreases. This is because of exhaustion of chemicals at the end of the run. This is a general trend seen in all runs, therefore you can remove these low quality bases from the end of your run. If you have incidentally a bad quality base in the middle of a read this is not the same general trend but is sporadic. Which is hard to automatically remove with tools such as cutadapt. If you really want an N for each low quality bases within the read, you can make a simple program for that, but why would you need that? Most aligners can handle a few bad quality bases within a read, and for variant calling you'll need many reads (so the incidental error base will be overruled by the good ones).


awk working with large files I have two very large vcf files 2GB and 6GB I want to look for unique combinations of CHROM and POS and output the row that matches. However, because the files a so large my machine always hang and stop processing. Is there are a way to work around the problem of this large file?

I am using this command, taken from an answer to another question:

awk  '{ if(NR==FNR){a[$1$2]=$0}else{if($1$2 in a){print}}}' file1.vcf file2.vcf

file 1

#CHROM  POS ID  REF ALT QUAL
1   10366   rs58108140  G   A
1   10611   rs189107123 C   G
1   51954   rs185832753 G   C
1   13327   rs144762171 G   C

file 2

#CHROM  POS ID  REF ALT QUAL
1   10366   rs58108140  G   A
1   51935   rs181754315 C   T
1   51954   rs185832753 G   C
1   52058   rs62637813  G   C
1   52144   rs190291950 T   A

output

1 10366 rs58108140  G   A
1 51954  rs185832753 G   C

Maybe comm, which is part of GNU coreutils, is more efficient.

comm writes to standard output lines that are common, and lines that are unique, to two input files; a file name of ‘-’ means standard input. Synopsis:
comm [option]… file1 file2

In your case you could sort both files based on the first 3 columns:
file1=file1.vcf
file2=file2.vcf

comm -12 <(awk '{print $1,$2,$3}' $file1 | sort) \
         <(awk '{print $1,$2,$3}' $file2 | sort) > matches

# result:
cat matches
1 10366 rs58108140
1 51954 rs185832753
#CHROM POS ID

This could be of course polished by removing the header.
Also if you want the original lines as output, you can grep the matches to the original files:
grep -f matches file1
#CHROM POS ID REF ALT QUAL
1 10366 rs58108140 G A
1 51954 rs185832753 G C

Edit: In the first version, intermediate temp files were used, but as suggested by terdon in the comments, this can be avoided using Process Substitution allowing the <(..) contstruct.


Filtering raw sequencing reads I have to filter the raw sequencing reads based on the following criteria:

Remove reads containing adapters
Remove reads containing N > 10% (N represents base that could not be determined)
Remove reads where the Qscore (Quality value) of over 50% bases of the read is <= 5

The reads are paired-end fastq files. I have been able to do the first part by using trimmomatic:

trimmomatic \
  PE in_1.fq.gz in_2.fq.gz \
  pair_out_1 unpair_out_1 pair_out_2 unpair_out_2 \
  ILLUMINACLIP:adapter.fa:2:30:5

The other 2 steps are a bit more challenging since I have not found a tool that perform exactly what is required. I could write some custom code of course but I would like to use the existing tools if possible. No need to reinvent the wheel.

Do you have any suggestions on how to proceed?

I think you can come close to this with cutadapt:

cutadapt \
  -a ${FWD_ADAPTER} -A ${REV_ADAPTER} \         ## adapter
  --max-n 0.1 \                                 ## do not allow > 10% Ns
  -q 5 \                                        ## remove bases with basequal < 5
  --minimum-length (readlength/2) \             ## minimum length after trimming
  -o read1_trimmed.fq -p read2_trimmed.fq \     ## outputs
  read1.fq read2.fq                             ## inputs

This will remove reads with more than 10% of bases being N and remove bases with base quality below 5. If you then set --minimum-length to half the read length, this should also remove reads where more than half of the bases had a base quality below 5 and therefore got trimmed away. Poor base quality in Illumina reads typically happens towards the 3' end so I think this little workaround is a sufficient proxy to achieve your goal. A base quality of 5 is very low in any case, so if sequencing went normal the typical quality should be notably higher. I would be surprised if this third filter is really necessary and will eliminate a notable number of reads.


Output file names are not known beforehand I have a folder where the hdf5 files inputs are (matrices of scRNAseq raw UMI counts). There is a python script that loads them, splits them into multiple matrices, normalizes them, and saves as csv files. When the name pattern of the output files is known, things are easy, I know how to do that, but when one input (hdf5 file) can produce various number of outputs (csv files) I am not sure how to approach the issue. Currently my Snakemake file relevant part looks the following way:

hdf5_files,=glob_wildcards('/main/{hdf5_file}.hdf5')

rule all:
  input:
--->  expand('data_files/normalized/{hdf5_file}.csv', hdf5_file=hdf5_files)
    ...

rule normalization:
  input:
     script = 'python/normalize.py',
     path = '/main/{hdf5_file}.hdf5'
  output:
--->  path = 'data_files/normalized/{hdf5_file}.csv'
  shell:
     "python {input.script} -data {input.path} -path {output.path}"

Above I highlighted the 2 lines that I am not sure now how I should write. Any suggestions would be greatly appreciated.

Could you give an example of names of input files and output files of normalize.py ?

It seems to me that your case could be similar to the one described about the dynamic output:

https://snakemake.readthedocs.io/en/stable/snakefiles/rules.html#dynamic-files

Here you would probably use it for the output of normalization and for the input of all (with a wildcard different from hdf5_file).

Disclaimer: I haven't had the opportunity to use this feature so far.


How to take into account alternative bwa mem mapping when computing coverage when mapping short reads with bwa mem, if a read has alternative mapping positions they are reported by bwa mem in the X0 and XA tags.
Now, let's say I want to compute the coverage of my bam file. How would I do for the alternative hits to be taken into accounts. What I mean is that, if a read QUERY has a selected best HIT position and 3 alternate HIT positions, I would like that QUERY read to be counted as 1 unit of coverage for the best HIT position, but also as 1 unit of coverage for each of the alternate HIT position.

My understanding is that samtools depth will only consider the best hit (which makes sense in 99.99% of the cases).

I guess I could edit the bam file with custom scripts so that the alternative hits are also reported as primary hits but I was wondering if there was an already existing tool (as I am not that confident with my ability to mess around with sam files).

Thanks a lot.

You would want to instead use the -a option with bwa mem to output all alignments. You would be best off filtering these with a custom script to only keep those with a close enough score to the primary hit. samtools depth can't be told to include secondary alignments (it's hard-coded), so you'll need to use something else like bamCoverage or samtools mpileup.


How to remove all BAM read groups from all reads (not just the header)? I have problem with one my BAMs---it appears to have invalid read groups.

Normally when I have such a problem, I remove all the read groups from the BAM header as follows:

samtools view -H your.bam | grep -v "^@RG" | samtools reheader - original.bam > new.bam

and then I would execute AddOrReplaceReadGroups from Picard:

http://broadinstitute.github.io/picard/command-line-overview.html#AddOrReplaceReadGroups

I have tried this, and there must still be reads with a read group in the optional tags RG:Z:.

How to do I definitively remove all read groups from a BAM, including all of the reads?

You might have to manually strip those auxiliary tags off:

samtools view -h your.bam | grep -v "^@RG" | sed "s/\tRG:Z:[^\t]*//" | samtools view -bo your_fixed.bam -

The sed bit is searching for the aux tag and removing everything up to the next tab.


STAR-long parameters for aligning RNA ONT reads to genome Are there any suggested parameters to align ONT reads to the reference genome using STAR-long? For now, I used the parameters suggested here, but I noticed a weird behaviour.

I have RNA reads (D. melanogaster) from R7 and R9 flowcells, separately. I only selected to analyze 2D reads in pass category.

I have, respectively, 113249 reads for R7 and 40318 reads for R9. I aligned those reads and get (only!) 150 uniquely mapped reads for R7 data and 8017 uniquely mapped reads for R9 data. I tried to run again the same command on a different server with fresh compilation, but the output file is consistent with these 150 reads.

However, if I align the same with GMAP, I get 78016 uniquely mapped reads for R7 and 33523 uniquely mapped reads for R9, so I suspect that something went wrong in the alignment run.

I am aware that the two mappers behave very differently, STAR-long being  more precise and preferring to report mappings of fewer reads but at better loci, and GMAP being overall less precise, trying to map the most of the reads but at not-so-good loci.

I was wondering if some of you had experience with this and could suggest me the best parameters for RNA reads from ONT?

I've had great results using minimap2, particularly when combined with a pre-treatment of Canu for error correction (using minimap2 for the read-to-read mapping):

# correct reads
~/install/canu/canu-1.6/Linux-amd64/bin/canu overlapper=minimap \
   genomeSize=100M minReadLength=100 minOverlapLength=30 -correct \
   -p 4T1_BC06 -d 4T1_BC06 \
   -nanopore-raw workspace/pass/barcode06/fastq_runid_*.fastq
# align reads
~/install/minimap2/minimap2 -p 10 \
   -a ~/db/fasta/mmus/ucsc/mmus_ucsc_all_cdna.idx \
   -x splice <(pv all/4T1_BC06/4T1_BC06.correctedReads.fasta.gz) | \
  samtools sort > 4T1_BC06_all_vs_mmusAll.bam

Update: the most recent nanopore basecaller, combined with the most recent minimap2, seems to now do a decent job with mapping, and pre-correction of reads no longer seems necessary. More recently I've been using LAST to map to the transcriptome, and minimap2 to map to the genome. Minimap2 has the ability to use a homopolymer-compressed genome index, which means that the most common consistent error for nanopore reads (i.e. misjudgement of homopolymer length) will not influence the mapping rate.


How to make Venn diagram for the Macs peak calling output of two data sets? I have two outputs of macs2 peak calling for two of my data sets. I wanted to plot the Venn diagram to see how many peaks are shared. I mainly work on Unix/Linux. Do you know any way that I can have those diagrams plotted?

Much Appreciated

Check out the tool Intervene. It's exactly designed for this. Github link.


Getting a stretch of genomic ranges from a dataframe/granges object based on metadata column I have a "test" data.frame object in R, which is basically a small subset of a 66000 row dataframe, which looks as follows:

chr   start    end      PC1     comp
chr1  3360001  3400000  -0.009  B
chr1  3400001  3440000  -0.004  B
chr1  3440001  3480000  -0.001  B
chr1  3480001  3520000  0.003  A
chr1  3520001  3560000  0.002  A
chr1  3560001  3600000  -0.001  B
chr1  3600001  3640000  0.002  A
chr1  3640001  3680000  0.004  A
chr1  3720001  3760000  0.003  A
chr1  3760001  3800000  0.003  A
chr1  3800001  3840000  0.007  A
chr1  3840001  3880000  0.006  A
chr1  3880001  3920000  -0.004  B
chr1  3920001  3960000  -0.017  B

As you can see, the interval between each start and end position is 40kb. The column comp can have only two values, A or B. Going from row 1 to the last row of my dataframe and based on the values of the column comp, I want to get the whole stretch/block of the genomic region, once the compartment(character in column comp) is changed, and perform a mean of the PC1 values in that stretch.

So in this example, I basically want to shrink my test dataframe as follows:

chr1  3360001  3480000  mean(c(-0.009,-0.004,-0.001))  B
chr1  3480001  3560000  mean(c(0.003,0.002))  A
chr1  3560001  3600000  mean(c(-0.001))  B  <---- #notice this is a 'singleton'
chr1  3600001  3880000  mean(c(0.002,0.004,0.003,0.003,0.007,0.006))  A
chr1  3880001  3960000  mean(c(-0.004,-0.017))  B

I tried with a simple code, but it only runs until the first compartment value change(so row 4 in this test), and I am unable to figure out how to make it run for the whole dataset of 66000 rows. Also the use of multiple for loops, which is what I tried, may not be the prettiest way to do this.

I guess there exists a Bioconductor package which performs this task very easily. If anyone can point me to that direction, it'd be great.

The GenomicRanges package in bioconductor will do this.  It lets you do 'set' operations on genomic ranges, such as concatenating them, finding overlaps, etc.  Hopefully something in the following code will be inspiring, even if I haven't solved the exact problem for you here:

library(GenomicRanges)

# make a GenomicRanges set corresponding to the original table
test.gr = makeGRangesFromDataFrame(df = test, start.field = "start", end.field = "end", seqnames.field = "chr")

# make separate GenomicRanges sets for each of the two partitions of the genome
test_A.gr = makeGRangesFromDataFrame(df = test[test$comp=="A",], start.field = "start", end.field = "end", seqnames.field = "chr")
test_B.gr = makeGRangesFromDataFrame(df = test[test$comp=="B",], start.field = "start", end.field = "end", seqnames.field = "chr")

# apply the reduce() function to concatenate segments which are adjoining each other into single segments
reduced_test_A.gr = reduce(test_A.gr)
reduced_test_B.gr = reduce(test_B.gr)

# group PC1 values by segment classification
test[queryHits(findOverlaps(test.gr, reduced_test_A.gr)), "PC1"]
test[queryHits(findOverlaps(test.gr, reduced_test_B.gr)), "PC1"]


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

Answer I received on Biostars that appears to work as expected:

bcftools norm -Ob --check-ref wx -f <reference_fasta> <vcf.gz> -o <output_name>


Tools to create annotated table of variants from VCF The problem: I have a VCF file, a reference genome, and a bunch of annotations for the reference (genes, repeat regions, etc.) as GFF or BED files.

What I would like is a tool that takes all of this as input and outputs a tab- or comma-delimited table containing as much information as possible. Potential columns in the output include:

Type of variant (SNV, structural, etc)
Details of variant (e.g. reference base, variant base, coverage, position, etc)
Annotations overlapping variant
Annotations near variant (e.g. is it just upstream of a gene)
If it appears in a coding region, does it change the amino acid.

There are many tools that do something like this. But to the newbie (like me) it is not clear which tools are worth starting with. Since the majority of tools take a little effort to get working in the first place, my question is:

What tool comes closest to doing what I have described, is known to reliable, and is likely to be maintained for the next few years.

Hail might be an option for you.

It is actively developed by a growing team at the Broad.  It is rigorously tested (continuous integration, continuous deployment, bug reports get regression tests, blah blah blah).

It was designed to solve this problem (among others).  It can import a variety of formats, including VCF, TSV, UCSC BED, JSON and interval files.  (We don't have explicit support for GFF, but we can probably handle them with general facilities.  If not, get in touch and we'll add support.)  It can call out to VEP (and soon Nirvana, Illumina's VEP rewrite).  It has general facilities to transform, filter, clean and query data.

What's more, we've curated a large collection of annotation resources (currently ~22 databases of annotations + VEP), hosted in a public bucket in Google cloud, and built an interactive query builder to select which resources you want to use.  Get in touch if you'd like us to add additional resources.  You can read more about it here.

Disclaimer: I work on Hail.


How exactly is "effective length" used in FPKM calculated? According to this famous blog post, the effective transcript length is:

$\tilde{l}_i = l_i - \mu$

where $l_i$ is the length of transcript and $\mu$ is the average fragment length. However, typically fragment length is about 300bp. What if when the transcript $l_i$ is smaller than 300? How do you compute the effective length in this case?

A related question: when computing the FPKM of a gene, how to choose a transcript? Do we choose a "canonical" transcript (how?) or combine the signals from all transcripts to a gene-level FPKM?

The effective length is $\tilde{l}_i = l_i - \mu + 1$ (note the R code at the bottom of Harold's blog post), which in the case of $\mu < l_i$ should be 1. Ideally, you'd use the mean fragment length mapped to the particular feature, rather than a global $\mu$, but that's a lot more work for likely 0 benefit.

Regarding choosing a particular transcript, ideally one would use a method like salmon or kallisto (or RSEM if you have time to kill). Otherwise, your options are (A) choose the major isoform (if it's known in your tissue and condition) or (B) use a "union gene model" (sum the non-redundant exon lengths) or (C) take the median transcript length. None of those three options make much of a difference if you're comparing between samples, though they're all inferior to a salmon/kallisto/etc. metric.

Why are salmon et al. better methods? They don't use arbitrary metrics that will be the same across samples to determine the feature length. Instead, they use expectation maximization (or similarish, since at least salmon doesn't actually use EM) to quantify individual isoform usage. The effective gene length in a sample is then the average of the transcript lengths after weighting for their relative expression (yes, one should remove $\mu$ in there). This can then vary between samples, which is quite useful if you have isoform switching between samples/groups in such a way that methods A-C above would miss (think of cases where the switch is to a smaller transcript with higher coverage over it...resulting in the coverage/length in methods A-C to be tamped down).


Can I change the name of file features.tsv to genes.tsv it is said that the features.tsv from Cell Ranger v3 is analogous to the genes.tsv from Cell Ranger v2. So can I change the file name to genes.tsv for Seurat to read it? I know that Seurat has the update version 3 that can read the features.tsv, and I would like to use the old version. So if I change the file name, did it cause any problem?

You can follow the 10X instructions to load your matrix into R:
the following is a direct quote from their webpage:

library(Matrix)
matrix_dir = "/opt/sample345/outs/filtered_feature_bc_matrix/"
barcode.path <- paste0(matrix_dir, "barcodes.tsv.gz")
features.path <- paste0(matrix_dir, "features.tsv.gz")
matrix.path <- paste0(matrix_dir, "matrix.mtx.gz")
mat <- readMM(file = matrix.path)
feature.names = read.delim(features.path,
                           header = FALSE,
                           stringsAsFactors = FALSE)
barcode.names = read.delim(barcode.path,
                           header = FALSE,
                           stringsAsFactors = FALSE)
colnames(mat) = barcode.names$V1
rownames(mat) = feature.names$V1

You can then create a seurat object using mat

object <- CreateSeuratObject(mat, project="")


How to map gene names to ensembl Ids in R I am trying to convert a list of gene names into Ensembl IDs but I'm having a bit of trouble.

I've tried using the EnsDb.Hsapiens.v86 package but I always get LRG-ID values returned as well, even when using the GeneIdFilter "startsWith" thing.
My code looks like this:

targets <- as.character(targetCSV[, 2])

targetEnsemble <- ensembldb::select(
      EnsDb.Hsapiens.v86,
      filter = GeneIdFilter("ENS", "startsWith"),
      keys = targets,
      keytype = "SYMBOL",
      columns = c("SYMBOL","GENEID"))

And for example I get back:

    SYMBOL          GENEID
1     CYR61 ENSG00000142871
2      CTGF ENSG00000118523
3    AMOTL2 ENSG00000114019
4    ANKRD1 ENSG00000148677
5    ANKRD1         LRG_379
6    IGFBP3 ENSG00000146674
7        F3 ENSG00000117525

I also tried using org.Hs.eg.db package but get a different problem where some of my genes return NA.

My code:

  mapIds(x = org.Hs.eg.db,
     keys = targets,
     "ENSEMBL",
     "SYMBOL",
     fuzzy = TRUE,
     multiVals = first)

And this works on some but not all the genes for some reason.

    $CYR61
[1] NA

$CTGF
[1] NA

$AMOTL2
[1] "ENSG00000114019"

$ANKRD1
[1] "ENSG00000148677"

$IGFBP3
[1] "ENSG00000146674"

$F3
[1] "ENSG00000117525"

If anyone can solve either of these problems that would be great, thanks in advance!

You can pull a full list of Ensembl gene IDs and MGI symbols from biomaRt in R. From there it should be easy enough to filter for the genes you want:

library(biomaRt)
hsapiens_genes <- getBM(attributes = c("ensembl_gene_id",
                                       "hgnc_symbol"),
                        mart = useMart("ensembl", dataset = "hsapiens_gene_ensembl"))

You do not find some genes in your example because you were not using the official gene name. For example CTGF is offocially called CCN2 and with this name you will find the appropriate Ensembl gene ID in the biomaRt table. This is why I personally prefer to use gene name that are like ENSG00000118523|CCN2 in my work. That way you always have a unique identifier (the Ensembl ID) that is unambiguous and you do not run into trouble that for one gene many more or less official gene names exist.


cn.MOPS fails with 'missing value where TRUE/FALSE needed' I am trying to use cn.MOPS to call CNVs on a set of whole exome sequencing data bam files. My script (up to the point where it fails) is:

source("https://bioconductor.org/biocLite.R")
library(cn.mops)
BAMFiles <- c('sample1.bam', 'sample2.bam', 'sample3.bam', 'sample4.bam', 'sample5.bam', 'sample6.bam', 'sample7.bam')
segments <- read.table("region.bed",sep="\t",as.is=TRUE)
gr <- GRanges(segments[,1],IRanges(segments[,2],segments[,3]))
X <- getSegmentReadCountsFromBAM(BAMFiles,GR=gr,parallel=11)
resCNMOPS <- exomecn.mops(X)

The last command, resCNMOPS <- exomecn.mops(X) fails with:

Error in if (all(segMedianT == 0)) { :
  missing value where TRUE/FALSE needed
In addition: Warning message:
In normalizeChromosomes(X, chr = chr, normType = normType, qu = normQu,  :
  Normalization for reference sequence  chrUn_gl000228 not applicable, because of low number of segments

If I run traceback immediately afterwards, I get:

> traceback()
4: FUN(newX[, i], ...)
3: apply(sINI[chrIdx, , drop = FALSE], 2, segment, minSeg = minWidth,
       ...)
2: cn.mops(input = input, I = I, classes = classes, priorImpact = priorImpact,
       cyc = cyc, parallel = parallel, normType = normType, normQu = normQu,
       norm = norm, sizeFactor = sizeFactor, quSizeFactor = quSizeFactor,
       lowerThreshold = lowerThreshold, upperThreshold = upperThreshold,
       minWidth = minWidth, segAlgorithm = segAlgorithm, minReadCount = minReadCount,
       useMedian = useMedian, returnPosterior = returnPosterior,
       ...)
1: exomecn.mops(X)

I am using R version 3.4.4 on an Ubuntu 16.04.5 LTS machine.

My bed file only has one region which is on chrUn_gl000228:

$ grep chrUn_gl000228 region.bed
chrUn_gl000228  112604  113879

If I remove that region and rerun the same code, it works as expected. So, my question is, how can I detect this and fix it automatically? Given a collection of bam files and the bed file describing the capture regions used to sequence the samples and generate the bam files, how can I check if cn.MOPS will complain about any of the segments?

Presumably, once I have built these two objects:

segments <- read.table("region.bed",sep="\t",as.is=TRUE)
gr <- GRanges(segments[,1],IRanges(segments[,2],segments[,3]))
X <- getSegmentReadCountsFromBAM(BAMFiles,GR=gr,parallel=11)

I should be able to then check them and see if they have the right data and, if not, remove the offending segment from the GRanges object. But how can I do this?

I managed to track down what was producing this error. The problem was that my bed file only had one segment for chromosome chrUn_gl000228. However, the cn.MOPS method relies on normalizing read counts across multiple segments for each chromosome. Since this chromosome only had a single segment, the exomecn.mops(X) was failing.

So, either remove the single segment from the bed file, or use a bed file that has at least 2 segments on each chromosome. I created a fake segment on chrUn_gl000228 and confirmed that the program worked as expected and didn't throw any errors.


Global alignment between two sequence X and Y with maximum number of identical matches If 2 protein sequence X of length m and Y of length n and if there are several highest scoring global alignments, then I want to get one alignment that has largest number of columns in which a letter from X is aligned to an identical letter from Y. I need an algorithm to find alignment with maximum number of identical matches. Can anyone help me?

There are a lot of possible ways to count the number of perfect matches in an alignment, and the details will depend on which alignment tool you are using.
But, if you use the built-in pairwise aligner in BioPython (Bio.Align.PairwiseAligner(), https://biopython.org/docs/1.75/api/Bio.Align.html), then the easiest way is to use the formatted string representation of the alignment, which uses the pipe/vertical-bar character ('|') to indicate perfect matches.
Example code:
from Bio import Align

aligner = Align.PairwiseAligner()
alignments = aligner.align("GAACT", "GAT")
max_matches = 0
for i, a in enumerate(alignments):
    align_str = f'{a}'
    align_str = align_str.split('\n')
    match_str = align_str[1]
    match_count = match_str.count('|')
    mismatch_count = match_str.count('.')
    print(f'Number of matches: {match_count}')
    print(f'Number of mismatches: {mismatch_count}')
    if match_count>max_matches:
        max_matches = match_count
        max_ind = i

best_alignment = alignments[max_ind]

In this example, both alignments have the same number of perfect patches (3), but that may not always be the case, depending on the input sequences and the scoring parameters used. The example just uses the default scoring parameters, but in general, you probably want to choose parameters appropriate for your type of alignment.


How to I get the infix expression for the body in a SBML function defintion using python libsbml? How do I translate the body of a function definition in SBML into infix using python libsbml?

For example after loading the model using:

reader = libsbml.SBMLReader()
doc = reader.readSBMLFromString (r.getSBML())
model = doc.model

To do the translation you need to access the function definition. First, make sure there are function definitions in your model by using:

if model.getNumFunctionDefinitions() > 0:
   # continue

To get the infix of the body for the ith function definition use the following code:

f = model.getFunctionDeinfition(i)
body = f.getBody()
print (libsbml.formulaToL3String(body))

As a side note you can also get the jth argument for the function definition using:

f = model.getFunctionDeinfition(i)
arg = f.getArgument(j)
print (libsbml.formulaToL3String(arg))


Exporting nt's differences from MEGA alignment I have an alignment of some full genome sequences of DENV1 that we have recently sequenced from the serum of patients. I have already cleaned the sequences, built the alignments and done the SNPs calling. After that I got the consensus sequences from each sample and I used MEGA to build the alignments of the consensus sequences and to construct a tree etc.

From just looking the alignment of the consensus sequences I can see several nt s differences among the sequences.

Does anyone know if you can export from MEGA like a sort of table where you could see the nt changes for each position and sample?

You don't say which version of MEGA you are using, but this works for MEGA X on Linux:

Open your alignment with the "Sequence data explorer"
At the "Highlight" menu, select "Variable sites".
Then at the menu "Statistics", select "Use only highlighted sites" and select what you consider an appropriate output - I choose "Display results as CSV".
Now at the menu "Data" select "Export data"; at the "Sequence data export options window" if necessary select the data format again, and select "For each site" at "Writing site numbers", and select "Only highlighted sites" at "Selected sites to include".

This will result in a comma-separated values file, showing only positions with variable sites.


Fitting branch models with godon Is it possible to fit branch models (Yang and Nielsen, 1998) in godon? That is, the model which averages over all sites, as opposed to the branch site models?

Branch-models are currently not implemented in godon. In the PAML manual Ziheng Yang writes:

This model is very parameter-rich and its use is discouraged.

This was one of the reasons I did not implement it.
Disclaimer: I am the author of godon.


Motif Discovery in Bipartite networks using Cytoscape I am trying to find over-expressed motifs in a bipartite transcription factor gene network.My genes are not very well characterized, so finding functional motifs has been a challenge. I used the motif-discovery app on Cytoscape and calculated the total number of 3 node motifs on the network. In my network, where transcription factors regulate genes, the only kind of three node motif that is possible is a simple input module.

This is an example of a simple input module where a transcription factor regulates multiple genes

I further created random networks to my real world network and statistically (Student's t-test) proved that there are more number of single input modules in the real rather than the random network. However, I am just not feeling satisfied with just this one explanation because I would have like to draw some biological inference from this. I was thinking of also analysis the 4 node and 5 node motifs in the network, find the number of times any of my hub nodes occupy a single node of these motifs and maybe use that and prove the importance of my hub nodes (though hub nodes are worthy of more exploration just by virtue of their high degree). Has somebody done something like this previously. It is a given that hub nodes would be part of the maximum number of motifs just by virtue of their high degree and connection to multiple genes. So would this prove to be redundant?

Thank you in advance for any help that could be provided.

I have recently develop an R package (sorry not a Cytoscape solution) to deal with this kind of situations of something grouping other elements, like pathways and genes.
The package is based on GSEABase of Bioconductor to store this relationships. Basically you should store a gene set for each TF and genes it is linked. I would recommend to use the name of the TF as the name of the gene set and then create a GeneSetCollection to be able to use the package.

The package provides with functionalities that allow you to calculate the distribution of genes per gene set (using genesPerPathway), as well as, pathwaysPerGene to know in how many other TF regulate each gene.

I am not sure of the biological inference you can take from this (I designed the package to test this) but you can also see if there is some specific preference for the TF to link to genes which are also regulated by other TF. This could be done by looking at the conditional probability (condPerPathways or sizesPerPathway).
Sorry the manual pages don't have examples and the vignettes are incomplete yet.

Also you could use the package tidygraph to calculate the centrality of the TF or of the genes you want.


Protein expression units? I'm starting for protein expression analysis.

In gene expression, we have TPM, RPKM etc. However, all I can see for protein expression is: "High", "Medium", "Low" and "Not detected".

  https://www.proteinatlas.org/ENSG00000140326-CDAN1/tissue

Q: Is there a unit for measuring protein expression like gene expression? I've seen people talking about PE (protein existence). Can we use it for protein expression?

I would guess that the real unit for measuring protein expression would be pM (pico-moles) or a count of protein molecules.

I'm not sure we have any methods to calculate such quantities though.


Customizing bigWig file I generate bigWig files using bamCoverage from deeptools, in part for my colleagues to visualize their mapped libraries in the IGV viewer.

A problem is that the displayed track name is apparently the file name, which is not convenient because some files for different libraries have the same name. The reason is that the files are organized with one directory per library. I would prefer to avoid changing the way I name and organize the files.

Is there a way to edit my bigWig files in order to have control of the track name independently from the file name?

As far as I understand, bigWig format derives from wiggle format, which have a track definition line that can contain a name and a definition attribute: https://genome.ucsc.edu/goldenpath/help/customTrack.html#TRACK.

Maybe there is a way to convert a bigWig file back to wiggle, edit the track definition, and then rebuild an updated bigWig file.

There's no equivalent to the wiggle header in bigWig (or bigBed) files, which is why UCSC uses the file name. This is actually the reason for the track line stuff that you linked to, since you can then specify a name and just point to where the bigWig (or other format) file is on the internet.

BTW, you can certainly convert your bigWig to wiggle, add the track line, and convert it back with UCSC tools (e.g., bigWigToWig and wigToBigWig). It shouldn't change the track name when you display it in UCSC (though if it does, then let me know, since I'll modify pyBigWig and libBigWig to support that then).


Bedtools wrongly indicates zero coverage I have a BAM file produced using BWA-MEM and I need to calculate the coverage of each exon in this sequencing. I tried using this command:

bedtools coverage -a exome_sorted.bed -b exome_alignment.bam -sorted

And the result was the following:

chr1    12595   12802   ref|DDX11L1,ref|NR_046018,ens|ENST00000518655,ens|ENST00000450305,ens|ENST00000456328,ens|ENST00000515242  .       ref|DDX11L1     0       0       207     0.0000000
chr1    13163   13658   ref|DDX11L1,ref|NR_046018,ens|ENST00000518655,ens|ENST00000450305,ens|ENST00000456328,ens|ENST00000515242  .       ref|DDX11L1     0       0       495     0.0000000
chr1    14620   15015   ref|WASH7P,ref|NR_024540,ens|ENST00000538476,ens|ENST00000438504,ens|ENST00000488147,ens|ENST00000541675,ens|ENST00000423562       .       ref|WASH7P      0       0       395     0.0000000
...

That is, zero coverage for all exons. But there are reads mapping to the exons, as I can see by opening my BAM file on IGV and looking at some exonic regions:

I have already tried sorting the BED file using sort -k1,1 -k2,2n (as recommended on the Bedtools documentation) and using bedtools sort with the -faidx option, but neither option solved my issue.

I have also checked if the chromosome names are the same in the BAM file and in the BED file, using samtools view -H. Both use chr1,chr2,...,chr22,chrX,chrY, although BAM has chrM and stuff like chrUn_gl000214 which are not on the BED.

Thanks!

EDIT:

@ATpoint pointed out that all visible reads in this region have MAPQ0. I looked at another exon that did not have this issue, but Bedtools still outputs zero coverage. IGV picture below:

If this is your output :

chr1    12595   12802   ref|DDX11L1,ref|NR_046018,ens|ENST00000518655,ens|ENST00000450305,ens|ENST00000456328,ens|ENST00000515242  .       ref|DDX11L1     0       0       207     0.0000000

Then I guess your bedfile looks like this (the first 6 fields) :

chrom : chr1
FeatureStart : 12595
FeatureEnd : 12802
name : ref|DDX11L1,ref|NR_046018,ens|ENST00000518655,ens|ENST00000450305,ens|ENST00000456328,ens|ENST00000515242
score : .
strand : ref|DDX11L1

The issue I see here is that you are giving the wrong information to the strand field of your bedfile. It should be either +, either - or . if you don't want to consider the strand information.

Hope it helps !


How to index residues by icode in biopython? I am indexing a structure by the code below:
PDB_CODE = '1lve'
RES_CHAIN = 'A'
RES_NUM = 27

structure = PDBParser(is_pqr=True).get_structure(PDB_CODE, PDB_CODE+".pqr")
residue = structure[0][RES_CHAIN][RES_NUM]

However, the 1lve contains residues with the same RES_NUM but different icodes 'A' to 'F' so my code either ignores them or throws an error if I pass '27A'. Does biopython allow this kind of structure indexing directly?

Okay, I found the answer myself. It is possible to index not only by RES_NUM, but by its id (check with residue.get_id()), so the below example would work:
#structure[0][RES_CHAIN][het, resseq, icode]
structure[0][RES_CHAIN][' ', 27, 'A']