What is a sensible number of gene/observations to explain PCA variance? I am working with a set of RNASeq dataset. I have about 4000 observations (genes) on 20 samples and plotting a PCA I found the clustering doesn't vary much when I use different number of genes, but the % variance of PC1 and PC2 does. I have ~1300 DEG. But how do you select a sensible cutoff for the number of genes/observations to explain PCA variance? Would it be correct to select the number of genes that are DEG which account for >50% of the variance between PC1 and 2?
Few examples:

All genes: PC1(34%), PC2(14%)

Top1000: PC1(45%), PC2(17%)

Top500: PC1(49%), PC2(19%)

Top50: PC1(55%), PC2(24%)

Top 5: PC1(75%), PC2(24%)

The general goal of PCA in RNA-seq can be stated as, "I'd like a low-dimension representation of my data to allow easy assessment of the gross structure of my samples, specifically for assessing missing batch effects, samples swaps, etc." Ideally, this low-dimension representation can fit into a plot or two. In other words, we'd like the first 2 (maybe 3) PCs to explain a reasonable amount of variation. A few things to note from this alone:

The number of differentially expressed genes were not mentioned. They have no relevance here.
There was no mention made of the treatment groups clustering together. Seeing this is NOT a goal of PCA.

As you have noticed, the more input rows (genes, assuming you haven't already transposed it) the less variation explained by the first 2 principal components (PCs). On the flip side, the more rows you input the better PCA is describing the structure of your entire dataset.
Question: Do I really care about the gross structure of the entire dataset?
Answer: No, I just want to see problems that I might need to account for.
So you really don't have to use all of the data in the PCA, just the bits that will allow you to see any issues. In general that'll be the few hundred most-variable genes (after some reasonable transformation so you're not just looking at the most highly expressed genes!). The exact number shouldn't really matter very much, something in the range of 300-500 most variable genes should suffice. A common practice is to just randomly pick a nice round number somewhere in or around this range. This is not a number that needs optimization unless perhaps you have a very large number of samples and therefore PC1 and PC2 only account for a small (say 10-20%) percentage of the variation, in which case you're better off looking at either a number of difference principal components or using a different dimension reduction technique, such as UMAP.


How can I do an overlapping sequence count in Biopython? Biopython's .count() methods, like Python's str.count(), perform a non-overlapping count, how can I do an overlapping one?

For example, these code snippets return 2, but I want the answer 3:

>>> from Bio.Seq import Seq
>>> Seq('AAAA').count('AA')
2
>>> 'AAAA'.count('AA')
2

For Biopython 1.70, there is a new Seq.count_overlap() method, which includes optional start and end arguments:

>>> from Bio.Seq import Seq
>>> Seq('AAAA').count_overlap('AA')
3
>>> Seq('AAAA').count_overlap('AA', 1, 4)
2

This method is also implemented for the MutableSeq and UnknownSeq classes:

>>> from Bio.Seq import MutableSeq, UnknownSeq
>>> MutableSeq('AAAA').count_overlap('AA')
3
>>> UnknownSeq(4, character='A').count_overlap('AA')
3

Disclaimer: I co-contributed the .count_overlap() methods with Peter Cock, see 97709cc


Analyzing transcripts for non-coding RNA From a RNA seq experiment we look for the counts of transcripts of tRNA but we are also interested of other types of RNA we can find. Because of the protocol used for the library prep we can't analyse miRNA but we can analyse ncRNA.
I never analyzed transcripts for ncRNA and I have some questions:

The name of the transcripts are like ENST00000611544.1, ENST00000387069.1, since they are non coding, are they linked to a gene? Should I rename the name of each ncRNA with a gene name or working with these long names is ok?
When analyzing transcripts of RNA I did enrichment analysis to see what pathways are more active. But when analyzing transcripts for ncRNA is possible to do something similar to enrichment analysis? How can analyse the functions of the ncRNA that we get from the study?

Thank you very much in advance

ncRNAs are by itself not special. Sure they might be more unstable and less abundant than other transcripts but for the analysis it is more or less the same. You should take your transcript level counts and aggregate them to the gene level. The tximport package comes to mind. Every transcript (given it is annotated) is linked to a gene. ENST00000611544 is connected to the U2 gene but these snRNAs are also short and I doubt you see them in normal RNA-seq.
In the end you can process data as any RNA-seq. Quantify your data against a transcriptome, e.g. with salmon, then aggregate transcript level abundance estimates (so basically counts) to the gene level and see whether your genes of interest have counts you can work with or whether you mostly have zeros. A transcriptome can be obtaine e.g. from GENCODE. The available annotation files will tell you which genes these transcripts originate from. Check for the GTF files.
Towards enrichment analysis, well it all comes down to the annotations you have. If the ncRNAs are annotated in any pathway then you can perform standard analysis with tools such as goseq or gprofiler2. If they are not annotated, well then there is no information available.


Adding external annotation to the phylogeny tree I was looking to overlay annotation on my phylogeny made,which was solved here question asked where i was looking to compare two different phylogeny which i could do.

Now i want to label the tips into various groups such as if its human or canine etc.One of the possible solution i found is this one not sure how to implement it.

alignment  file
annotation file

library(dendextend)
library(seqinr)
library(phytools)
library(phangorn)

a<-read.alignment("clean_dup_align_fast.fas", format="fasta")
a.phydat<-as.phyDat(a)
dist.a.phydat<-dist.dna(as.DNAbin(a.phydat))
upgma.a<-upgma(dist.a.phydat)
parsimony(upgma.a,a.phydat)
pars.a <- optim.parsimony(upgma.a, a.phydat)
pars.a<-acctran(pars.a, a.phydat)
pars.a.rooted<-root(pars.a, outgroup="AAA64460", resolve.root=T)
pars.a.rooted.dd<-as.dendrogram(force.ultrametric(pars.a.rooted))

i want to label according to the start column which is my source in my annotation file
Any suggestion or help would be really appreciated

My working solution using this , this works as I did the alignment of sequences using full header name which kind of help to look for the pattern and grep it.
The data file

I would still like to see the ggtree way instead of putting lots of label in the plot I would only like to see if I can annotate all the tips which have one organism source as one color rather putting all the labels

pars.a
# "Phylogenetic tree with 516 tips and 514 internal nodes"

# by using the class() function
class(pars.a)
# "phylo"

# or by using the str() structure function
str(tree)
str(pars.a)
# "List of 4"
# this list includes $edge, $Nnode, $ tip.label and $edge.length
# the tree$tip.label includes family designation
#tree$tip.label  # 516 of these
pars.a$tip.label
# from the Science paper, we have seven kinase families:
# kinase categories... TK, TKL, STE, CK1, AGC, CAMK, CMGC
# with the following colours
# "red", "green", "paleblue", "orange", "yellow", "purple", "pink", "green"

# by using the grep()function on the tree$tip.label part of the object
# we can find the tip labels that include "TK/" - i.e. tyrosine kinases
#grep("TK/", tree$tip.label)  # gives a list of numbers with "TK/" in tip label
#length(grep("TK/", tree$tip.label))

grep("Human",pars.a$tip.label)
length(grep("Human",ml.a.rooted$tip.label)
)
# thus there are 94 tip labels with that are designated TK (not TKL tyrosine kinase like)

# make a vector for each tip.label called tipcol with black on all of these...
#tipcol <- rep('black', length(tree$tip.label))
tipcol <- rep('black',length(pars.a$tip.label))

# make a vector with our list of kinase categories
#kinaseCats <- c("TK/", "TKL", "STE", "CK1", "AGC", "CAMK", "CMGC", "RGC")
kinaseCats <- c("Human","Rat","Murine","Bovine","Canine","Feline","Porcine","Gallus","Bat","Infectious","Transmissible")

# make a vector of color we want:
colorsList <-c("red", "darkolivegreen3", "blue", "orange", "blueviolet", "purple", "pink", "green","brown","navyblue","purple4")

# replace colours where grep gives "TK" as red, etc in a loop
#for(i in 1:length(kinaseCats)){
 # tipcol[grep(kinaseCats[i], tree$tip.label)] <- colorsList[i]
#}

for(i in 1:length(kinaseCats)){
  tipcol[grep(kinaseCats[i], pars.a$tip.label)] <- colorsList[i]
}

# plot with edge length false to see nodes better
#plot(tree,
 #    use.edge.length = FALSE,
  #   tip.color=tipcol,
   #  cex = 0.25)

plot(pars.a,
     use.edge.length = FALSE,
     tip.color=tipcol,
     cex = 0.5)
# slow to draw due to text - a bit annoying!
nodelabels(cex=0.4)

figure

there are nice features in R ggtree package. You may use the geom_cladelabel() function for your purposes. But I haven't find yet a way to automatize clade annotation. To do it, you have to be sure of the clade homogeneity before removing tiplabs. here are some commands to do clade label annotation:

annot<-read.delim("annota.txt", header = T)

ggtr<-as.treedata(pars.a)
g<-ggtree(ggtr)
#find the node number corresponding to the common ancestor of the top clade in you tree
MRCA(ggtr, "AAL40400", "AAF25519")
#here, check that all accessions belong to the same organism/group
subset(annot, label=="AAL40400" | label=="AAF25519")
#in this case, sequences are related to bovine.
g<-g + geom_cladelabel(node=89, label="Bovine", align=F, fontsize=4)
plot(g)

#another example
MRCA(ggtr, "AAA46456", "AAF97738")
subset(annot, label=="AAA46456" | label=="AAF97738")
g<-g+ geom_cladelabel(node=73, label="Rodent", align=F, fontsize=4)
plot(g)


Intersecting a data frame based on the multiple columns I have 2 data frames

> head(df1)
   Chromosome Start_Position End_Position       Gene
1:      chr12       46254706     46254706       ARID2
2:      chr17        7579849      7579849        TP53
3:       chr9       21974777     21974780      CDKN2A
4:       chr1       27105930     27105930      ARID1A
5:       chr1      201982973    201982973        ELF3
6:       chr3       89480459     89480459       EPHA3
> head(df2)
  Chromosome Start_Position End_Position
1          1          10583       863511
2          1       12841835     12854479
3          1       17231299     17232877
4          1      120380739    142540462
5          1      142540484    142684781
6          1      142685823    143462598
>

How I can intersect genes from df1 which are in the range of df2? Because both data frames have chromosome number, start and end but I don't know how to identify genes from df1 based on df2

I have tried

df1[df1$Start_Position%in%df2$Start_Position,]

But I obtain one gene

I also tried

> new_dataset <- df1 %>% right_join(df2, by=c("Chromosome","Start_Position"))
Error: Can't join on 'Start_Position' x 'Start_Position' because of incompatible types (character / integer)
>

reformat df2 so that you have chr1 is instead of 1, then you basically have two bed files.  Then use bedtools intersect.  https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html


How to use the DI-tector software for defective genomes analysis? I am trying to use a free .py script from this site

  http://www.di-tector.cyame.eu/

while unfortunately there's neither manual nor any documentation presented.

I have downloaded the .py script however I have absolutely no idea how it can be started. Can anyone give me an idea of how that can be done?

The corresponding article (https://pubmed.ncbi.nlm.nih.gov/30012569/) only displays parameters for executing it and there is no way to understand what steps must be performed firstly. Thank you.

See the help:

python DI-tector_06.py -h
usage: DI-tector_06.py [-h] [-g HOST_REF] [-s MIN_SEGMENT] [-m MIN_MAPQ]
                       [-n NB_READS] [-o OUTPUT_DIRECTORY] [-t TAG] [-d]
                       [-l INDEL_LENGTH] [-f] [-p POLARITY] [-q] [-k]
                       [-x NB_THREADS]
                       Virus_Ref Input_Data

positional arguments:
  Virus_Ref             Virus genome reference sequence in FASTA format.
  Input_Data            File containing single reads in FASTQ format.

optional arguments:
  -h, --help            show this help message and exit
  -g HOST_REF, --Host_Ref HOST_REF
                        Host genome reference sequence in FASTA format.
  -s MIN_SEGMENT, --Min_Segment MIN_SEGMENT
                        Minimum segment length. Default is 15.
  -m MIN_MAPQ, --Min_MAPQ MIN_MAPQ
                        Skip alignments with MAPQ smaller than INT. Default is
                        25.
  -n NB_READS, --Nb_Reads NB_READS
                        Show only DVGs with counts reads > or egal to INT.
                        Default is 2.))
  -o OUTPUT_DIRECTORY, --Output_Directory OUTPUT_DIRECTORY
                        Enter a directory name that all compiled output files
                        will be saved in.
  -t TAG, --Tag TAG     Enter a tag name that will be appended before each
                        output file. Default is 'DI-tector'.
  -d, --DVG_sequences   Generate multi-fasta file with DVG sequences. Default
                        is (OFF).
  -l INDEL_LENGTH, --InDel_Length INDEL_LENGTH
                        Skip alignments with size of InDels smaller or egal to
                        INT. Default size is 1.
  -f, --Fasta           Select '-f' if data is in FASTA format fasta. Default
                        is FASTQ.
  -p POLARITY, --Polarity POLARITY
                        [0] Positive strand genome / [1] Negative strand
                        genome. Default is 0.
  -q, --No_Quantification
                        Inactive percentage quantification. Quantification
                        need bedtools Default is (ON).
  -k, --Keep_files      Keep intermediaite files (i.e. alignment etc...).
                        Default is (OFF).
  -x NB_THREADS, --Nb_threads NB_THREADS
                        Number of threads. Default is 1.


What are the state-of-the-art cell-type RNA-Seq deconvolution methods? I would like to find the proportion of each cell-type in bulk RNA-Seq transcriptomics data. I am looking for some guidance on the following:

What are the state-of-the-art methods?
What are their assumptions? (e.g. do they assume that cell-type expression signatures are available in advance?)

I see a lot of people using xCell. There should also be papers systematically compare many different methods, which might interest you.


Ribosomal RNA QC quantification using SILVA I would like to create a simple QC check of RNA-Seq data that simply maps the data to rRNA and then counts the number of reads that map.

I've done this manually for human using sequences for 5S, 5.8S, 18S, and 28S rRNA (rDNA, actually).  I gathered those sequences manually to create a small fasta file to use as a reference and it worked well.

I would like to have a more definitive source for this reference though, and I don't necessarily want to use a species-specific database, so I thought that using a fasta file downloaded from SILVA would be great, but I don't know if I'm using the wrong fasta file or what.  I'm not very familiar with SILVA, but I noticed that when I tested it using the same sequencing data, I got unexpectedly many fewer hits than I did with my manually created 4-sequence reference than with the SILVA NR99 fasta file.

So I blasted my 4 rDNA sequences against the SILVA NR99 data and only the 18S was found.  The others (5S, 5.8S, and 28S) don't appear to be there.  Am I misunderstanding the nature of the SILVA project or could I be doing something wrong?

UPDATE: Please note, this question is about where to find rRNA and/or rDNA sequences to use as a reference, not about tools that can be employed in quantifying rRNA read content.  While I appreciate hearing about tools that can generate nice reports about rRNA reads in an RNA-Seq run, the question is not answered unless it addresses the missing rRNA/rDNA sequences that are necessary for comprehensive quantification.  Any tool can count rRNA/rDNA mappings, but the results will only be as good as the rRNA reference data supplied (unless it is a reference-less tool).  And it's the lack of 5S, 18S, and 28S ribosomal subunits in SILVA that is my main concern.  It's the difference between 1k reads and 500k reads mapping to rRNA/rDNA sequences.

SILVA is intended for metagenomic samples and thus does not contain the human ribosomal subunits.  While the 5.8S subunit is in the ensembl, ncbi, and ucsc annotations, and there are numerous rRNA "biotype" (or "biomol") entries, the 5S, 18S, and 28S subunit sequences are absent.  Even downloading rRNA entries in the rmsk table from ucsc will not return all 4 major human subunits.

You can however (as I discovered) download all rDNA gene sequences (including those 4 subunits) from the nucleotide db on ncbi using the search term "biomol_rrna[PROP]".  This returns 240,240 coding sequences (or 30,775 gene sequences).  You can then use this as a reference for mapping, blasting, or as input to rRNA QC tools.

The following e-Utilities commands from ncbi will download the gene sequences:

esearch -db nuccore -query "biomol_rRNA[prop]" | efetch -format fasta > all_rrnas.fasta

This is not a common analysis.  Usually, rRNA reads can just be ignored and analyses are unaffected, which is likely why there is a paucity of records for these sequences in annotation files, but mapping to just these sequences could make for a simple QC measure as a gauge for rRNA contamination.  While the sequences may be present in the genome, they're not in the annotations and may even be masked as repetitive, since there are so many copies.  So you have to be careful to count unique reads and not mappings.


Trimmomatic trims most of the reverse strand I have a paired-end data. It is from small part of human genome not a WGS or WES. I use BWA for alignment and do variant calling. Since I had some false positives I wanted to do trimming.
I used trimmomatic with command:
trimmomatic PE -phred33 -threads 2 read1.fastq.gz read2.fastq.gz read1_trimmed.fastq read2_trimmed.fastq LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15`

I do not know why but I have this output:
Input Read Pairs: 1834513
Both Surviving: 1822034 (99.32%)
Forward Only Surviving: 5740 (0.31%)
Reverse Only Surviving: 417 (0.02%)
Dropped: 6322 (0.34%)
TrimmomaticPE: Completed successfully

I do not have any adapters since the are removed in bcl2fastq step. In FASTQ quality control I have seen that quality of reverse and forward reads are similar in quality but trimmomatic trimmed many reads in reverse file. I cannot do alignment with these files, BWA gives a warning that says:
[bseq_read] the 2nd file has fewer sequences.
[process] read 155012 sequences (16693324 bp)...
[bseq_read] the 2nd file has fewer sequences.

Is there anyone who faced a problem like this before?
Thank you in advance

Maybe you need to remove the unpaired reads:
trimmomatic PE -phred33 -threads 2 \
    read1.fastq.gz read2.fastq.gz \
    read1_trimmed.fastq.gz read1_unique.fastq.gz \
    read2_trimmed.fastq.gz read2_unique.fastq.gz \
    LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15

Then you can try with bwa.


scRNA-seq, 10x cellranger pipelines I'm starting to do sc-rnaseq using 10x cellranger pipelines, and i add TdTomato sequence to mouse reference genome and add an entry in the gtf. My code is But when I makref, it remindered that:
Writing genes GTF file into reference folder...
Invalid number of columns in GTF line
Can anyone tell me the reason?

As was suggested by @benn, you will need to add the TdTomato sequence to the end of your fasta formatted genome file. If your TdTomato sequence is in fasta format already, cat will do the trick. If not, just go the end of the file, start a new line and add a header line starting with ">" then immediately the name of your gene, no whitespcaes in between, something like ">TdTomato". Then add your sequence in a new line.

The format of the GTF file is explained in the image you posted, the trick is to use the exact same name (in column 1) that is used in the genome file you will create at the previous step (no ">" though). Again, this line, in which the different fields (columns) will be separated by tabs, will be added to the end of your GTF file.


creating graph of distribution of blast hits on the query sequence I'm attempting to visualize the results of my BLAST search in a way similar to the graphical display of the distribution of blast hits from the web BLAST.

For example, from my BLAST search:

except I would like the graph to contain information including information for the distribution of all 6139 hits across my query sequence.

Since this is a large number of sequences to visualize a distribution, I am considering using a score of the number of hits on a certain region, and having an output more similar to:

If I were to run my blast using the BLAST+ command line tool, and read the results into R, which parts of the output should I plot in order to recreate the graphical display of the distribution? Ideally I would like to make a reusable object that can make an equivalent plot for any BLAST input.

Thanks.

You probably want to include query start (qstart) and query end (qend) in your blast output.

Something like this:

blastn -query your.fasta -out blast.out.txt -db your.db -outfmt '6 qseqid sseqid qstart qend length evalue'

In R you can take the "qstart:qend" from each line for density plot.

There are many ways in R to plot the densities of these start and end amino acids.

Let me show an example with a small data frame:

qstart <- c(200, 300, 250, 400, 500)
qend <- c(300, 450, 400, 600, 650)

df <- as.data.frame(cbind(qstart,qend))

aa <- vector()
i=1
for(i in 1:5){
aa <- append(aa, c(df[i,1]:df[i,2]))
i+1
}

hist(aa)
dens <- density(aa)
plot(dens)


How to download sequence alignments against UniProtKB sequence database for a particular Pfam Id using R? I am trying to download sequence alignments for families against the UniProtKB sequence database as shown in the example below using R code. Can someone kindly suggest me how to?
https://pfam.xfam.org/family/PF02171/alignment/uniprot/format?format=pfam&alnType=uniprot&order=t&case=l&gaps=default&download=1

The following code served my purpose.
library(readr)
pfam_url <-
  'https://pfam.xfam.org/family/PF02171/alignment/uniprot/format?format=pfam&alnType=uniprot&order=t&case=l&gaps=default&download=1'
seqs <-
  read_table(
    pfam_url,
    col_names = FALSE,
    col_types = cols(X1 = col_character(), X2 = col_character())
  )


Delete all 4 lines of a fastq read from a fastq file using read ID I have the following error when running bowtie2:

  Error: Read HWI-D00466:116:CC62WANXX:3:1102:7363:63646 1:N:0:GCACACG
  has more read characters than quality values.

I now want to remove all 4 lines of this specific read from the fastq file.

How can I use awk or sed to do this?

Don't do it: Your FASTQ file is either malformed or your FASTQ record spans more than four lines, which is allowed in FASTQ.  For a detailed description of what can go wrong in FASTQ parsing see for example http://biopython.org/DIST/docs/api/Bio.SeqIO.QualityIO-module.html#FastqGeneralIterator.

If the FASTQ is malformed, then you should really ask yourself how this happened in the first place and fix the source of the problem. If the record is valid FASTQ, then I suggest parsing the read with for example FastqGeneralIterator and dumping the parsed result back to FASTQ in a 4-line-per-record form.


How do you calculate the top L/5 score? I am currently investigating how to predict protein structure and contact map predictions. There, I see things like top L/5 and top L/2 scores as a way to evaluate the contact map. What formulas do you use to evaluate these? Also, what is a meaningful evaluation method?

A contact map predictor will usually output the probability of a contact for every residue pair in a protein, often ommitting very close pairs (e.g. less than 6 residue separation).
To calculate the L/n score, rank the predicted contacts descending by probability and look at the first L/n contacts. The fraction of these that are actually contacts in the native structure is the L/n contact precision score.
Often, you would filter the list to look at only contacts of a certain sequence separation, e.g. long range contacts between residues of >23 separation (definitions vary between papers). If you have a good predictor this score is high for L/10 and falls as you go to L/1, which is expected as proteins don't usually have that many long range contacts to find.


Extracting expression data from GSE dataset downloaded from GEO I have downloaded GSE16146 dataset from GEO using GEOquery R package.

I would like to extract "Data table" from downloaded GSE16146.

>library("GEOquery")
>GSE16146 <- getGEO("GSE16146")
>Table(GSE16146)[1:5,]

This returns the following error:

> Table(GSE16146)[1:5,]
Error in (function (classes, fdef, mtable)  :
unable to find an inherited method for function ‘Table’ for signature ‘"list"’

Can anyone help?

According to the manual, all you need to do is:

library('GEOquery')
gseGSE16146 <- getGEO('GSE16146', GSEMatrix=FALSE)

As explanation, getGEO() outputs by default to GSEMatrix=TRUE and returns a list of ExpressionSet objects. You should get what you were looking for with:

Table(GSMList(gseGSE16146)[[1]])[1:5,]

The manual has also a paragraph about this:

  The GSE has  a  metadata section, just like the other classes.
  However, it doesn’t have a GEODataTable.  Instead, it contains two
  lists, accessible using GPLList and GSMList, that are each lists of
  GPL and GSM objects.

I tested the code I posted on R 3.3.2 and I get access to the Data Table just fine.

update:

Answering how to extract the actual expression data (not what asked in the original question):

library(GEOquery)
data = getGEO("GSE16146")
datExpr = exprs (data[[1]])


Y Chromosome Aligned Reads in scATAC-seq data from a female-derived cell line? I'm working with scATAC-Seq data on the K562 cell line, which is supposed to be derived from a female patient. While following the scATAC-seq data analysis pipeline, after performing bowtie alignment they recommend filtering out all reads aligned to mitochondrial DNA and on the Y chromosome.

Out of curiosity, I decided to count how many reads aligned to chrY and found that it can be quite high - chrY has approximately 10% as many reads as a similarly-sized chromosome (chr19 in hg19).

Is it normal to find reads mapping to chrY when the cell line is female-derived, or did I mess up somehow?

There are homologous regions between X an Y chromosomes: https://en.wikipedia.org/wiki/Pseudoautosomal_region

It is therefore normal to have some female-derived reads mapping in Y chromosome.

You should probably check what proportion of such reads fall in other parts of the Y chromosome than pseudoautosomal regions.


Write the output of the VariantAnnotation package to file I'm having some problems in exporting the annotations of some variants from a VCF file, after using the "VariantAnnotation" package from bioconductor. Below the script I've used:

#!/usr/bin/env Rscript

library(VariantAnnotation)
library(BSgenome.Scerevisiae.UCSC.sacCer3)
library(TxDb.Scerevisiae.UCSC.sacCer3.sgdGene)

args = commandArgs(trailingOnly=TRUE)

vcf <- readVcf(args[1], "sacCer3")

sc <- BSgenome.Scerevisiae.UCSC.sacCer3
txdb <- TxDb.Scerevisiae.UCSC.sacCer3.sgdGene

coding <- data.frame(predictCoding(vcf, txdb, seqSource=sc))

The coding dataframe looks like this:

  seqnames start end width strand paramRangeID REF     ALT QUAL FILTER
1     chrI   455 455     1      +         <NA>   C       T   40   PASS
2     chrI   455 455     1      +         <NA>   C           40   PASS
3     chrI   455 455     1      +         <NA>   C       G   40   PASS
4     chrI   455 455     1      +         <NA>   C       A   40   PASS
5     chrI   456 456     1      +         <NA>   T       A   40   PASS
6     chrI   456 456     1      +         <NA>   C CCTGCTC   40   PASS
  varAllele CDSLOC.start CDSLOC.end CDSLOC.width PROTEINLOC QUERYID TXID CDSID
1         T          121        121            1         41       1    1     1
2                    121        121            1         41       2    1     1
3         G          121        121            1         41       3    1     1
4         A          121        121            1         41       4    1     1
5         A          122        122            1         41       5    1     1
6   CCTGCTC          122        122            1         41       6    1     1
   GENEID    CONSEQUENCE REFCODON  VARCODON REFAA VARAA
1 YAL069W  nonsynonymous      CTC       TTC     L     F
2 YAL069W not translated      CTC        TC
3 YAL069W  nonsynonymous      CTC       GTC     L     V
4 YAL069W  nonsynonymous      CTC       ATC     L     I
5 YAL069W  nonsynonymous      CTC       CAC     L     H
6 YAL069W  nonsynonymous      CTC CCCTGCTCC     L   PCS

But if I try to save it to file, for instance using the following command:

write.table(coding, 'test.tsv', sep='\t', row.names=F)

The output table looks like this:

"seqnames"      "start" "end"   "width" "strand"        "paramRangeID"  "REF"   "ALT"   "QUAL"  "FILTER"        "varAllele"     "CDSLOC.start"  "CDSLOC.end"    "CDSLOC.width"  "PROTEINLOC"    "QUERYID"       "TXID"  "CDSID" "GENEID"        "CONSEQUENCE"   "REFCODON"      "VARCODON"      "REFAA" "VARAA"
"chrI"  455     455     1       "+"     NA      "C"     <S4 object of class "DNAStringSet">     40      "PASS"  "T"     121     121     1       41      1       "1"     1       "YAL069W"       "nonsynonymous" "CTC"   "TTC"   "L"     "F"
"chrI"  455     455     1       "+"     NA      "C"     <S4 object of class "DNAStringSet">     40      "PASS"  ""      121     121     1       41      2       "1"     1       "YAL069W"       "not translated" "CTC"   "TC"    ""      ""

Notice the <S4 object of class "DNAStringSet"> bit. Is there any way to force R to give me the same representation as head? I couldn't find anything in Biostrings' documentation.

Is there any way to force R to give me the same representation as head?

head just cuts off everything after the first six entries. The magic, as it were, happens in the default print method for the specific type. But the representation generated by it is only useful for ad-hoc display, not for storage in a structured data file.

The best solution is therefore to manually convert the column types. Since it’s only one column:

coding$ALT = sapply(coding$ALT, as.character)

I’m a bit puzzled that the REF and ALT column have different types, anyway: that’s probably a bug in the {VariantAnnotation} package that should be reported.

Furthermore, you probably want to use quote = FALSE for the write.table. In general, write.table has terrible defaults. readr::write_tsv is almost always a better choice.


How is a principal component analysis for population structure made from a bam or other alignment file? I'm curious how a PCA for population structure is built from an alignment/BAM or SAM file. How do programs like angsd and ngsCovar actually build a PCA? What exactly makes things come out as different and how can this be determined? If I have 3 alignments that look like this:

Locus1:

>Indiv1
ACGT
>Indiv2
ACGC
>Indiv3
ACGA

Locus2:

>Indiv1
GACT
>Indiv2
GACA
>Indiv3
GAGT

Locus3:

>Indiv1
GGTC
>Indiv2
GGCT
>Indiv3
GCTA

(These example alignments leave out quality scores for simplicity, assume that each has high quality across the board.)

What algorithm or system would be used to graph them on a PCA? How many dimensions would that PCA have? Is it possible to model this by hand so that I can understand what's going on/how it is built?

In general, the process works as such:

Variants are found at a number of loci in each sample (missing loci are often imputed).
Each sample is scored at each locus, often according to the number of copies of the reference allele that it has (so, 0-2).
The resulting matrix is used for PCA.


How to search for high coverage SRA entries The Question
I want to find high coverage SRA entries, e.g., above 100x.
I guess the best way is to use https://www.ncbi.nlm.nih.gov/sra with an appropriate search term. I don't mind if the search results contain some "false-positives". (I.e., if some of the found entries have coverage above 80x, while I aimed for 100x, that's fine.)
I want to avoid using Entrez.esummary on a very large number of SRA entries (e.g., through Biopython), as it would be significantly slower (if I understand correctly).
How can I do that?

My Attempt
I tried using the Mbases field.Assuming the reference genome is at most 100kbp (e.g., in case I am interested in viruses), I thought a good approximation would be to search for SRA entries with more than 100*100k bases, i.e., at least 10M bases.
According to https://www.ncbi.nlm.nih.gov/books/NBK3837/, I thought I could use the search term 10:9999999[Mbases], but disappointingly, this didn't work. Just to make sure my syntax is right, I verified that 100:101[ReadLength] works. So it seems that Mbases doesn't support the range syntax.

If processing a huge table with lots of "false positives" is OK, then you could enter your search term, from top right choose "Send to" > "File" Format: "RunInfo".


What is a site in a phylogenetic tree? When I read about a term related to phylogenetics, for example, $K_a$/$K_s$ ratio, I constantly see the word site, whose meaning isn't explained and it puzzles me.

For example substitutions per synonymous site, substitutions per non-synonymous site, substitutions at a site etc.

Could you define a phylogenetic site?

So there were 0.6 substitutions per the blue-violet site, what does it mean? What is the length of the site? Is the length of sites constant? Is it the length of a sequence read? Are introns the part of the site?

A site is simply an individual discrete position, normally a single nucleotide (or amino acid). For example, consider a toy alignment:

sequence1 ATGC
sequence2 ATCC
#           ^ arrow points to the third site

Here there is 1 mutation in 4 sites, so 0.25 substitutions per nucleotide site.

I'm not sure where the term originated, but probably one of the earlier definitions is from renowned population geneticist Motoo Kimura, who simply stated in 1968:

  In this paper, “site” refers to a single nucleotide pair

A pair in this context, means either a homozygous or heterozygous site.

Synonymous sites are positions where nucleotide changes lead to no amino acid change, these are typically at the 3rd position of the codon where the nucleotide is often degenerate, for example, TAT, TAC TAG, and TAA all encode the Serine amino acid. By contrast non-synonymous sites are where a single nucleotide change does lead to an amino acid change, and therefore possibly a more important functional change to the resulting protein.

Thus, the ratio of Ka/Ks (or Dn/Ds) is thus often used as a crude proxy of selection, with Ka/Ks < 0, Ka/Ks = 0 and Ka/Ks > 0 implying purifying (negative) selection, neutral evolution and positive selection respectively. Note that this assumes the neutrality of synonymous sites, which is increasingly questioned due codon usage bias and regulatory elements overlapping coding regions (and also ignores mutational biases and biased gene conversion, although these can be accounted for in an evolutionary model, and also assumes no sequencing errors and perfect alignments, which can impact a lot on the results too).


Estimating BAM file from compressed fastq file size Is there a way to estimate the size of a BAM file will have after mapping with BWA?

The input file are two mates fastq files, compressed with gzip, each one about 70G.

A BAM file should have a roughly similar size to a compressed FASTQ file, because they're both compressed files that contain similar information. BAM will be slightly larger, because it contains information about mapping and differences from the reference.

If you are concerned about space and have a reference sequence available, you can store mapped reads as CRAM format, which will save a lot more space.


Generate VCF from different .bam files with different chromosome names I have two resources of .bam files. One is generated by our lab (1 sample = 1 bam). One is downloaded online (again 1 sample = 1 bam). For the downloaded samples the chromosomes are labelled: chr1, chr2, chr3 etc For our lab samples, the chromosomes are labelled: 1, 2, 3 etc. I want to generate a single VCF file of variants across all samples.

I'm using bcftools:

bcftools mpileup -Ov -f ref.fasta -b samples.txt | bcftools call -mv -o bamMge.vcf

However, I get no calls and the repeated error: [E::faidx_adjust_position] The sequence "1" was not found

I have two questions:

Is my strategy correct (i.e. bcftools mpileup)? Do I need to incorporate extra steps (NB I also have matching gVCF files)
How can I either (i) alter the chromosome labeling of one of the subsets of bam files, or, (ii) use a mapping file to match chromosome labels during the mpileup run?

Thanks

The main issue here is that your bam files have different chromosome labels. (i.e. 1,2,3 vs. chr1,chr2,chr3) as you mentioned.

This suggests that the data was mapped to a slightly different version of the reference genome. As a first layer, you might want to double check that these reference chromosomes are cross compatible (e.g. both correspond to hg19 for instance). To find more information on the reference genome used, you can check the header information using "samtools view -H "

Assuming that the reference genomes are cross compatible, you can try converting it to sam, adding the missing "chr" label in a pipe, and converting it back to a bam file. However, this will also require modification of the header file which can be challenging in itself.

The easiest (and more rigorous) approach is probably just to convert the batch of bam files with the wrong chromosome labels back into fastq files. You will then want to remap them using the exact same reference genome as the other batch of bam files. After which, you should be able to perform variant calling on the combined batch of bam files.


What are Orthogonal DNA Sequences? When reading a paper, I came across the following statement

  Standard vectors were constructed, flanked by pairs of homology sequences derived from yeast barcodes at the ends of each expression cassette. We reasoned that as these barcode sequences were designed to be orthogonal, they could serve a dual purpose of reducing the probability of mis-annealing and dictating the assembly order of multiple cassettes.

I am not sure what the "orthogonal" means here. I think it is about synthetic biology.

Orthogonal DNA sequences are the shortest possible single stranded nucleotide stretches constructed. These sequences are used as probes in DNA hybridization. Because of its length and construction it has least chance of binding with other probes.

This technique improves the efficacy and specificity and is clinically used to detect genetic disorders in shortest time duration.

Via:
https://www.nature.com/articles/srep24650

https://www.pnas.org/content/106/7/2289


Find a map/correspondence between two versions of a genome I am working with two versions of the C. elegans genome. I am finding interesting regions (specifically, tRNA genes) in version 1 and then I would like to know if version 2 also has a tRNA gene in that region. For this, I need to know the correspondence between coordinates in genome 1 and coordinates in genome 2. I think a whole-genome alignment programme might help me, e. g. nucmer or mauve. I have run nucmer on the two genomes and got a delta file but I don't know how to interpret it. I have run show-coords on the delta file and obtained a .coords file. It seems more human-friendly. The first few lines look like this:

/Users/user/Documents/project/output/genome/celegans/ref/wormbase/seq/c_elegans.PRJNA13758.WS245.genomic.fa /Users/user/Documents/project/output/genome/celegans/strain/seq/bristol/scaffolds.fa
NUCMER

    [S1]     [E1]  |     [S2]     [E2]  |  [LEN 1]  [LEN 2]  |  [% IDY]  |  [LEN R]  [LEN Q]  |  [COV R]  [COV Q]  | [TAGS]
===============================================================================================================================
       1      432  |      498      929  |      432      432  |   100.00  | 15072434 15371751  |     0.00     0.00  | I  I
       1      432  |     1002     1433  |      432      432  |   100.00  | 15072434 15371751  |     0.00     0.00  | I  I

Does this mean that chromosome I starting at 1 and ending at 432 of version 1 corresponds to chromosome I starting at 498 and ending at 929 and that there is 100% identity between these two 432-long sequences? In this case, it's easy to find a correspondence between bases in version 1 and bases in version 2. However, some lines do not align sequences of the same length:

   1   237139  |     4022   241181  |   237139   237160  |    99.99  | 15072434 15371751  |     1.57     1.54  | I  I

Also, in this sequence position 1 of chromosome I is aligned again? So is it aligning to more than one position in version 2 of the genome? Also one sequence is 237139 bases long and the other 237160 bases long. How do I know the internal structure?

Can anybody guide me, please?

I've also aligned the two genomes using Mauve:

progressiveMauve c_elegans.PRJNA13758.WS245.genomic.fa scaffolds.fa --output scaffolds.mauve > scaffolds.txt

I've tried to import the .mauve file into Python using Biopython but I got an error:

alignment = SeqIO.read('scaffolds.mauve', "mauve")

Traceback (most recent call last):
  File "code/mauve.py", line 4, in <module>
    alignment = SeqIO.read('scaffolds.mauve', "mauve")
  File "/Users/user/anaconda/lib/python3.6/site-packages/Bio/SeqIO/__init__.py", line 712, in read
    raise ValueError("More than one record found in handle")
ValueError: More than one record found in handle

Note that the alignment opened in the Mauve GUI looks fine. Unfortunately, the picture won't load here.

I am working on Mac OSX.

I would appreciate your help in understanding the delta file, resolving the Python error or any other support in finding a function between coordinates in version 1 of the genome and coordinates in version 2 of the genome.

A relatively convenient method would be to create a .chain file, as used by tools like liftOver and crossmap (and probably remap from NCBI). For your use-case, you can simply use the DoSameSpeciesLiftOver.pl command as described in the link. The resulting chain files are pretty self-explanatory and can be used to map BED files of regions of interest.


Why are inversions defined as the reverse complement and not just the reverse of the reference? I can’t quite understand the way inversions are defined. In particular I expect an inversion to be only a reversed version of the reference and not its reverse complement. Most sources use diagrams such as below which totally obscure the fact that it is also complemented.

Could someone provide an explanation on what causes an inversion and why it is defined the way it is?

If an inversion were just reversed then there would be 3' -> 3' bonds and 5' -> 5' bonds. That's why inversions are reverse complemented, you then maintain the normal 5'->3' direction. This is nicely covered over on our biology sister site, which also includes this image:


Is there a way to import tables in PubMed Central to MySQL or other Relational databases? I need to use these PMC tables from the medical journals for the purpose of integrating the data with a decision-making software. Though, I can see that one can access full-text for an article from PMC here (which includes the article table), but is there a way to extract the table from this XML document and export it to a local MySQL database or any other relational database?

If you can use R, consider using web scraping.

See here : http://bradleyboehmke.github.io/2015/12/scraping-html-tables.html


Can pair end reads with high heterozygosity be used to polish PacBio assembly? I have used canu (correct) and smartdenovo to assembly PacBio reads. Next I am going to polish my assembly using illumina pair end reads by Pilon. However, I found there is high heterozygosity in my pair end reads when I used jellyfish to do k-mer analysis. There are two peaks in the k-mer frequencies distribution plot. The size of my genome is about 190M, and it is a diploid organism.

What I am doubting is that can I use this pair end data to polish my assembly?

If not, I have to do sequencing again.

Resequencing your sample is not going to make it any less heterozygous. I am assuming that you are working with a diploid organism? Please update your question and specify the genome size and ploidy of your organism (the type of organism will help us update our answers, too).

The answer, yes, you can use Pilon + short reads to correct the PacBio assembly even with a highly heterozygous sample. To do so you should use the --diploid option in Pilon to let it know that there will be heterozygous sites. The way that it corrects heterozygous SNPs and heterozygous indels is different, and it seems to perform best with SNPs and sometimes does not correctly identify heterozygous indels (source).

Pilon itself doesn't handle phasing, but just corrects with whatever reads are in the bam file that it is provided with. The end result of providing a bam with reads from both haplotypes is that the corrected sequence will be a mosaic of the dominant allele at the heterozygous site as was observed in the short read data. This may or may not select for some alleles based on sequencing coverage bias of Illumina sequencers.


Books on application of Machine Learning in Bioinformatics I'm looking for books on applications of machine learning and statistical data mining in bioinformatics with example codes provided in R and/or python.
I came across the following two books: the first being a quite good one, provides the code in MATLAB, which I'm not very familiar with. The second one too is a good book with examples in R.

Feature Selection and Ensemble Methods for Bioinformatics by Oleg Okun

Data Science and Predictive Analytics: Biomedical and Health Applications using R by Ivo Dinov

Can people suggest other books based on their first hand experiences?
Thanks

One could read Bioinformatics: The machine learning approach as an overview of the established applications of machine learning to sequence alignment, RNA structure prediction, etc...
...Yet, this book is 20 years old, and there have been many important developments in bioinformatics and machine learning in this period, which have not yet been covered in any books. In other words, reading books on bioinformatics and/or machine learning is good for getting firm basic background in either of them, but any such books are inevitably out of date at the time of their publication.
I would like however to stress the importance of the basic background - one regularly sees announcements about breakthroughs achieved by reinventing techniques that were already described in the book that I mentioned.


How can I speed up INDEL calling/correction on BAM files? The samtools mpileup command has quite a neat feature that it is able to correct mapping errors associated with misalignment of INDELs. By default, the mpileup command will not work for reads that have more than 250X coverage of the reference genome. While this limit can be increased, very high coverage causes the mpileup program to grind to a halt, so it'd be nice to know if there's some easy way to make that faster.

To add a bit more context, I've been doing this with mitochondrial genome reads that were extracted from both Illumina whole-genome sequencing (coverage ~1000X), and from targeted amplicon sequencing done on the IonTorrent (coverage up to ~4000X).

I see that @rightskewed has mentioned the downsampling ability of samtools with samtools view -s <float> (see here), which seems like it might work as a solution for this if used prior to the mpileup operation.

I wasn't aware of the samtools subsampling when I had this problem a couple of years ago, so ended up writing my own digital normalisation method to deal with mapped reads. This method reduces the genome coverage, but preserves reads where coverage is low.

Because I was working with IonTorrent reads (which have variable length), I came up with the idea of selecting the longest read that mapped to each location in the genome (assuming such a read existed). This meant that the highly variable coverage for different samples (sometimes as low as 200X, sometimes as high as 4000X) was flattened out to a much more consistent coverage of about 100-200X. Here's the core of the Perl code that I wrote:

  if(($F[2] ne $seqName) || ($F[3] != $pos) || (length($bestSeq) <= length($F[9]))){
    if(length($bestSeq) == length($F[9])){
      ## reservoir sampling with a reservoir size of 1
      ## See https://en.wikipedia.org/wiki/Reservoir_sampling
      ## * with probability 1/i, keep the new item instead of the current item
      $seenCount++;
      if(!rand($seenCount)){
        ## i.e. if rand($seenCount) == 0, then continue with replacement
        next;
      }
    } else {
      $seenCount = 1;
    }
    if(($F[2] ne $seqName) || ($F[3] != $pos)){
      if($output eq "fastq"){
        printSeq($bestID, $bestSeq, $bestQual);
      } elsif($output eq "sam"){
        print($bestLine);
      }
    }
    $seqName = $F[2];
    $pos = $F[3];
    $bestLine = $line;
    $bestID = $F[0];
    $bestFlags = $F[1];
    $bestSeq = $F[9];
    $bestQual = $F[10];

The full code (which works as a filter on uncompressed SAM files) can be found here.


about dotplot legend meaning Here is the code I used:

DotPlot(object = cluster, genes.plot = c("Hcn4","Shox2","Nkx2-5","Nppa","Tnnt2","Gja5","Tbx18","Tbx3","Isl1","Scn5a","Cacna1d"),  plot.legend=TRUE)

I was confused about the legend with the Negative number shown for avg.exp.scale. Does it mean that the cluster expressed the genes at low level or that gene was down regulated? (like the Shox2 in the "0" compared with "2").

The expression values for each gene are scaled / standardized by subtracting the genes mean expression and dividing by its standard deviation. A value of -1 would imply it's one standard deviation below the mean expression for that gene.


Detect transcript isoform abundance for a specific gene in scRNA-seq I want to detect the count of isoform transcripts for a specific gene in scRNA-seq data. Data is coming from cells of Mus Musculus.
For transcript isoforms I mean the different alternatives provided by known annotation. For each one I need the total transcript count in each cell.

By looking at the count matrix, alternative transcripts are not defined (see for instance this dataset).

Should I come back to the fastq? Which tool may I use to produce a count matrix with counts per each alternative transcript?

Transcript quantification is a difficult enough problem as it is, when you add the extra difficulty of going from the low read numbers available in scRNAseq if gets even more difficult.

Added to this, many scRNAseq techniques (including 10X, which is what the linked dataset is) only capture the 3' end of the gene, so you can distinguish transcript with different 3' end, but the data simply doesn't exist for doing much more with standard tools.

Lior Patchner does claims that there is benefit in looking at "transcript compatibility counts" rather than just gene counts, and that you can use his Kallisto tool for this, and the these 3' end tagging tools do give you some data for the rest of the transcript.

You could also isolate transcripts with differing three 3' ends and map reads to these using the something like the UMI-Tools pipeline.

Alternatively you could look at things on a junction level - you won't capture all junctions, but you may be able to look at differential junction usage for those you do capture using something like SUPPA or MISO. For this I would use UMI-Tools to get de-duplicated BAM files, split them on cell identity and then pass them into your tool of choice. I don't really how much you will capture.

The last thing to say is that this might be a space worth watching for future developments (cough).

Data from some other single cell methods, like SMART-seq2 or C1 will sequence the whole transcript and you may be able to apply standard tools such as Salmon, Kalisto or RSEM to get transcript usage for them, although as I said, transcript level quantification is noisey, single cell data is noisey, single-cell level transcript data will be very noisy.

Whatever, you are going to have to go back to the raw sequencing data.

Disclosure: I am the author of UMI-Tools.


wtdbg2: practical implications of k-mer fsize and psize choice I am using wtdbg2 2.3 to assemble a human genome (sequenced on PromethION from a cell line). I filtered out reads with low average quality, and now I am trying to determine the parameters that will optimize overall assembly quality.

The wtdbg2 assembler has two k-mer related settings, documented as follows.

 -k <int>    Kmer fsize, 0 <= k <= 25, [0]
 -p <int>    Kmer psize, 0 <= p <= 25, [21]
             k + p <= 25, seed is <k-mer>+<p-homopolymer-compressed>

The requirement that k + p is <= 25** suggests that these parameters can be set simultaneously. However, in all of the presets for the various long read sequencing technologies, one of the two settings is always set to 0.

My questions are as follows.

Are k and p ever intended to be used simultaneously?
If so, how does the behavior change with -k 2 -p 19 versus -k 19 -p 2?
I'm familiar with some common implications of choosing k size for k-mer based analyses of Illumina reads. (Smaller k-sizes are more sensitive but also more subject to spurious matches. Larger k-sizes are more specific but also more likely to include a sequencing error.) Do similar implications apply for -k and -p with wtdbg2?
BONUS: Is there an explanation anywhere of how -k and -p are used? It looks to me like -k is related to k-mers in the traditional sense, and -p is a homopolymer compressed k-mer. Are seeds made by concatenating the regular k-mer and the compressed k-mer? It's possible I missed some important details in the preprint but I haven't been able to make the connection yet.

**Incidentally, the documentation here is wrong. During runtime wtdbg2 enforces a constraint of <= 23, not 25. This matches the description in the bioRxiv preprint stating that 46 bits (2 bits per base) are used to store each k-mer.

There is some description of homocompressed polymers (-p kmer) in this paper: https://arxiv.org/pdf/1708.01492.pdf


Convert Cotton Probe ID to Gene Symbol I am new to bioinformatics, my background is in Electrical Engineering. I am trying to convert Affymetrix Cotton Probe IDs to gene symbols. I have a gene expression dataset and I need the expressions for only certain genes. So in the dataset, I have a probe id such as Ghi.10301.1.S1_s_at, how can I get the gene symbol associated with it? Thanks

EDIT: This is the specific dataset I am looking at
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE50770

Your dataset was produced on a microarray. The probe IDs in the FASTA file are similar to read IDs produced from other NGS experiments like RNA-Seq. According to the workflow provided by Ensembl the first step is to align your reads to a reference. You'll only keep the reads that match the reference completely. Once the probes are aligned you'll be able to quantify the number of reads mapping to the gene that the probe was designed to attract.

I did find a file providing the database entries associated with each probe ID (one of which is the ID you provided) for the Affymetrix Cotton Genome Array.


Parallel version of cn.mops' bamDataRanges I am trying to run cn.MOPS to detect CNVs in multiple bam files. Following the example code, I have run:

source("https://bioconductor.org/biocLite.R")
biocLite("cn.mops")
options(width=75)
set.seed(0)
library(cn.mops)
library(Biobase)
library(GenomicRanges)
library(GenomeInfoDb)
cn.mopsVersion <- packageDescription("cn.mops")$Version
library(cn.mops)
BAMFiles <- c('bam1.bam', 'bam2.bam', ..., 'bamN.bam')
bamDataRanges <- getReadCountsFromBAM(BAMFiles)

PLEASE BE PATIENT... this might take a while. Consider using the parallel version of this function

Since this is indeed taking a while, I'd like to follow the tool's suggestion and use "the parallel version of this function". But how do I do that? I don't see any mention of "parallel" in the manual. Is there some standard Bioconductor/R trick to make a function parallel? Am I missing something in the manual?

How can I run bamDataRanges parallelized?

This is an option for the function documented in its help page. Just run ?bamDataRanges and you will see:

Usage:

     getReadCountsFromBAM(BAMFiles, sampleNames, refSeqNames, WL = 25000,
       parallel = 0, ...)

[...]

parallel: The number of parallel processes to be used for this
          function. Default=0.

So all you need to run with 10 cores is:

bamDataRanges <- getReadCountsFromBAM(BAMFiles, parallel=10)


How to interpret amino acid representation What is the difference between C383S and C383A? For C383S, does it mean Cysteine(C) position 383 in the sequence and its next amino acid connect to it is Serine(S)?

This means that at position 383 the Cysteine(C) is mutated to a Serine(S).  Depending on the organism and gene usually the mutation notation is linked to a specific transcript since the position could vary depending on transcript splicing.


How to make a distinction between the "classical" de Bruijn graph and the one described in NGS papers? In Computer Science a De Bruijn graph has (1) m^n vertices representing all possible  sequences of length n over m symbols, and (2) directed edges connecting nodes that differ by a shift of n-1 elements (the successor having the new element at the right).

However in Bioinformatics while condition (2) is preserved, what is called a De Bruijn graph doesn't seem to respect condition (1). In some cases the graph doesn't look anything like a de Bruijn graph at all (e.g. http://genome.cshlp.org/content/18/5/821.full).

So my question is, if I want to make it explicit that I am using the Bioinformatics interpretation of a de Bruijn graph, is there a term for it? Something like "simplified de Bruijn graph", "projection of a de Bruijn graph", or "graph of neighbouring k-mers"? Are there any papers making this distinction, or did I get it all wrong?

Several papers have made this distinction, and a few indeed use different terms to distinguish between them. For example, Kazaux et al. (2016) acknowledge that:

  These constraints favour the use of a version of the de Bruijn Graph (dBG) dedicated to genome assembly – a version which differs from the combinatorial structure invented by N.G. de Bruijn.

Kingsford et al. (2010) also recognise the distinction:

  Note that this definition of a de Bruijn graph differs from the traditional definition described in the mathematical literature in the 1940s that requires the graph to contain all length-k strings that can be formed from an alphabet (rather than just those strings present in the genome).

The oldest reference I found for a specific term to refer to the assembly-related structure is Skiena and Sundaram (1995), where they call it a subgraph of the de Bruijn digraph. Later, in 2002, Błażewicz et al. will refer to it as a de Bruijn induced subgraph.  The term de Bruijn subgraph is also formally defined in Quitzau’s thesis (2009). There, and also in the article (Quitzau and Stoye, 2008) the authors describe the sequence graph as a modification of the sparse de Bruijn subgraph (commonly used in assembly problems), where non-branching paths are replaced by a single vertex. The term sparse de Bruijn graph is also used by Chauve et al. (2013).

Another term that I found was word graph, described by both Malde et al. (2005) and by Heath and Pati (2007) as a subgraph or as a generalization of a de Bruijn graph. Rødland (2013) summarises some of the terms used for this data structure:

  The data structure is best understood in terms of the de Bruijn subgraph representation of S[k]. (...) Some authors may refer to this as a word graph, or even just a de Bruijn graph.

Although we can recognise that the distinction is not very relevant, the question is asking specifically for the situation where one wants to make such a distinction.


In calculating the retention index, why do we use the character state with the lowest frequency? To calculate the retention index for a phylogenetic tree, we use the following formula:

$$\frac{\text{maximum number of steps on tree - number of steps on the tree}}
{\text{maximum number of steps on tree - minimum number of steps in the data}}$$

To calculate the maximum number of steps, we look at the number of observed character states for a character and take the lowest number. For each character, we take the lowest number of observed traits and then sum across all characters.

Take the following dataset (columns are characters, taxa are rows):

          123456
Outgroup  000000
Taxa A    100001
Taxa B    110001
Taxa C    111001
Taxa D    111101
Taxa E    111111
Taxa F    111110

The maximum number of steps for the first character is 1, since the character state 0 is attested once and character state 1 occurs six times.

What I don't understand is why we do not use 6 for the maximum number of changes for the first character.

The maximum number of steps (or changes) is the number of taxa with state 1 or 0, whichever is smaller. Conceptually, this is the number of character changes if each taxon evolved its state independently of the other taxa.

If for a given character the character state 0 occurred twice and character state 1 occurred five times, 2 would be the maximum number of changes.


Difference between PPI networks and Gene Co-expression networks What exactly is the difference between protein-protein interaction (PPI) networks and gene co-expression networks? Based on my rudimentary understanding, PPIs are constructed from gene co-expression networks, and so represent processed data. But could someone elaborate more on the differences?

Based on my rudimentary understanding, PPIs are constructed from gene co-expression networks

Not really, no. PPI stands for protein-protein interaction. In PPI networks, each node is a protein and each edge (line) connecting it to another node represents an interaction. A very common example is that an edge represents a physical interaction between these two proteins. This means that there have been experiments to indicate that these two proteins come into contact in the cell and interact.

For example, they might form part of a protein complex, so their interaction was seen using chromatin immunoprecipitation, where the two proteins were captured together. Or, they may have been seen together using FRET, a technique that lets you detect when two molecules are in very close proximity.

Alternatively, some interactions can be theoretical, where in silico modelling has suggested that these two proteins can interact, even though we don't have any direct evidence that they do so.

Finally, yes, it is also possible to make a PPI network based on gene co-expression data, under the assumption that if two genes are co-expressed, then their protein products will be present in the same tissue. However, that doesn't mean they interact in any way, so it isn't really a very common way of building PPI networks, as far as I know.

In co-expression networks, on the other hand, each node is a gene and not its protein product. And here, the edges indicate the coexpression. Often the level of correlation between the two genes' expression data.

The most important difference between these two kinds of network, apart from the obvious mentioned above, is that they are used in different ways. A gene coexpression network lets you identify genes that are regulated in similar ways. It is a tool you can use to investigate the regulatory network of the cell, to see what genes are "turned on" or "off" together.

The PPI network, on the other hand, is a way of understanding the result of this gene co-regulation. You are not looking at how the gene network is regulated any more, but how the products of the genes interact with one another. So it is a very different focus.

You can maybe think of this in terms of school children and their parents. You could make a network where each node represents a pair of parents, and then draw lines between them if their children were born in the same year. That would give you a general overview of when this group of people procreated. However, that two children were born on the same day doesn't necessarily imply that the children themselves will ever come into contact. This could be the equivalent of the gene co-expression network.

Then, make another network where each node is a child and draw a line between any children who know each other. This would be the equivalent of the PPI network. You would likely find some correlations between the two, with children who know each other being more likely to have been born at around the same time, but the two networks give you very different information.


Can I submit a R package to Bioconductor or CRAN if I have already published it a journal? I have written a bioinformatics package in R that I want to publish in a bioinformatics Journal. Presently, I am maintaining a local repo of that package and I want to put in the Bioconductor repository (after publication or at least after submission to the journal), is it possible? or I must submit it to Bioconductor first?
Thanks

You can submit a package to a repo after the publication, but I would say that you really should do it before.

Bioconductor accept only packages that are not published on CRAN, however, academic publications are fine.

Conversely, Bioinformatics does not enforce any particular platform for sharing the code, but you have to make it available somehow (e.g. GitHub repo, university webpage, bioconductor). From the instructions for authors

  If the manuscript describes new software tools or the implementation
  of novel algorithms the software must be freely available to
  non-commercial users at the time of submission, and appropriate test
  data should be made available.

However, you as a developer should make all the effort to make your software available and easily accessible. Also, sooner you make it available, sooner people start to cite you, which is a great advantage for you. Finally, you really need users to make sure that your package really works (i.e. debugging all the corner cases).

I could not emphasise more, the correct order is 1. make it available 2. write a publication


SNP data Download Help I am new to the field of Bioinformatics. Single Nucleotide Polymorphism or SNP data are available on the the Internet to download.SNP related data are valuable for variety of reasons ranging from  disease related charatertcics and  precision medicine to name the few.

Can somebody explain me the steps of downloading this kind of data from NCBI or other resource.
For example if I want to download the the SNP in bunch of humans on a particular Gene from the specific start to end. What are the sequence of steps I would do online to download such data?

Update

I want SNP data for example for a particular gene say X from positions in range say base pairs 100,000 to 200,000 for individuals say 100 in txt or csv format or standard format that is available?

Like if I choose Gene X and position ranging from 100 to 109 base pairs in 5 persons I want following kind of data in txt or csv format with say two SNV at position 101 and 109

A A A A A G G G G G
A G A A A G G G G A
A A A A A G G G G A
A A A A A G G G G A
A G A A A G G G G G

regards

I would recommend to do this in R, using biomaRt.

For example if you want all the SNPs from IDH1 (ENSG00000138413) gene:

library(biomaRt)

ensembl<- useMart("ensembl")
ensembl <- useDataset("hsapiens_gene_ensembl",mart=ensembl)

IDH1 <- getBM(attributes=c('ensembl_gene_id','variation_name'), filters = 'ensembl_gene_id', values ="ENSG00000138413", mart = ensembl)

IDH1

You can add much more information, take a look at the help page of biomaRt.


Converting Gene Symbol to Ensembl ID in R I'm trying to convert ~20,000 different human gene symbols to ensembl IDs. I've been trying to use biomaRt to do this, but continue getting the following error

getBM( attributes=c("ensembl_gene_id") , filters= "mgi_symbol"  ,mart=ensembl)
Error in martCheck(mart) :
  No dataset selected, please select a dataset first.  You can see the available datasets by using the listDatasets function see ?listDatasets for more information.  Then you should create the Mart object by using the useMart function.  See ?useMart for more information

Is there a different format I should be using for this? I tried using BioConductor's description of how to use getBM but am a R novice and confused as to how to set this up to convert gene symbols into ensemble IDs.

You can use the standard annotation package for humans (you can also use biomaRt, but it can be more confusing, see below):

library("org.Hs.eg.db") # remember to install it if you don't have it already
symbols <- mapIds(org.Hs.eg.db, keys = ensemblsIDS, keytype = "ENSEMBL", column="SYMBOL")

The ensemblsIDS should be a character vector: c("ENSG001", "ENSG002", ..)

The solution with biomaRt:

library("biomaRt")
ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl")
getBM(attributes='hgnc_symbol',
      filters = 'ensembl_gene_id',
      values = ensemblsIDS,
      mart = ensembl)


Biopython Genbank writer not splitting long lines I am parsing a csv file of annotated sequences and using Biopython to generate Genbank files for each. I want to add annotations of the sequence features.

My output file shows features listed without the correct line breaks. Other software is then unable to parse the names of the features

misc_feature    58..60  /note="Oxidation"
misc_feature    100..102  /note="Oxidation"
misc_feature    241..243  /note="Oxidation"

It should be formatted like this

misc_feature    58..60
                /note="Oxidation"
misc_feature    100..102
                /note="Oxidation"
misc_feature    241..243
                /note="Oxidation"

My section of code adds a space between the location and the /note="" but the Biopython writer doesn't add line breaks.

for i in range(len(start_list)):
    feature = Feature()
# Strict Genbank format calls non-standard things misc_feature
    feature.key = "misc_feature"
    liab_type = "{}".format(liability).split(":")[1]
#need to add 1 to start list since csv counting is off
    feature.location = "{}..{}".format(int(start_list[i]) +1, end_list[i]) + "  " + '/note=\"{}\"'.format(liab_type)
    container.features.append(feature)
    i +=1

The Biopython source is toggling the wrapping of the feature qualifiers on no_space_key and Bio.GenBank._BaseGenBankConsumer.remove_space_keys

How do I enable wrapping of text to format output correctly?

    def __str__(self):
    """Return feature as a GenBank format string."""
    output = Record.INTERNAL_FEATURE_FORMAT % self.key
    output += _wrapped_genbank(self.location, Record.GB_FEATURE_INDENT,
                               split_char=',')
    for qualifier in self.qualifiers:
        output += " " * Record.GB_FEATURE_INDENT

        # determine whether we can wrap on spaces
        space_wrap = 1
        for no_space_key in \
                Bio.GenBank._BaseGenBankConsumer.remove_space_keys:
            if no_space_key in qualifier.key:
                space_wrap = 0

        output += _wrapped_genbank(qualifier.key + qualifier.value,
                                   Record.GB_FEATURE_INDENT, space_wrap)
    return output

I read the Biopython Genbank record source code carefully and realized I wasn’t using the Qualifier class correctly.

After reading the source code again, I realized I was using the data structure incorrectly.  I need to use the Qualifier class to add a key and feature.  This is formatted correctly and solves the problem.

feature = Feature()
# Strict Genbank format calls non-standard things misc_feature
feature.key = "misc_feature"
#need to add 1 to start list since csv counting is off
feature.location = "{}..{}".format(int(start_list[i]) +1, end_list[i])
liability_qualifier = Qualifier()
liability_qualifier.key = '/note="'
liability_qualifier.value = '\"{}\"'.format(liability).split(":")[1]
feature.qualifiers.append(liability_qualifier)


how to generate structure file from blast results of an input sequence I have tried to incorporate blast in my application using biopython.i could extract the blast results but now i want to extract the structure file corresponding matches. My code is :

from Bio.SeqIO import parse
from Bio.SeqRecord import SeqRecord
from Bio.Seq import Seq

file = open("prats.fasta")

records = parse(file,"fasta")
for record in records:
    print("Id: %s" % record.id)
    print("Name: %s" % record.name)
    print("Description: %s" % record.description)
    print("Annotations: %s" % record.annotations)
    print("Sequence Data: %s" % record.seq)
    print("Sequence Alphabet: %s" % record.seq.alphabet)

from Bio.Blast import NCBIWWW
sequence_data = open("prats.fasta").read()
sequence_data

result_handle = NCBIWWW.qblast("blastn", "nt", sequence_data)
result_handle

with open('results.xml', 'w') as save_file:
    blast_results = result_handle.read()
    print(blast_results)
    save_file.write(blast_results)

from Bio.Blast import NCBIXML
E_VALUE_THRESH = 1e-20
for record in NCBIXML.parse(open("results.xml")):
    if record.alignments:
        print("\n")
        print("query: %s" % record.query[:100])
        for align in record.alignments:
            for hsp in align.hsps:
                if hsp.expect < E_VALUE_THRESH:
                    print("match: %s " % align.title[:100])

from here i want to generate the xml files and then fetch structure of each nucleotides results

If you the DNA/RNA crystal structure:

Simply change the nt in NCBIWWW.qblast("blastn", "nt", sequence_data) to pdbnt to get the nucleotide sequence in a chain.
The PDB code will be the accession ID with an underscore plus chain letter.
Do note that the PDB structure may be offset.

If you want the crystal structures of the encoded protein:

blastx (translation) on the pdbaa database. As above the PDB code will be as above.

If you want a model of the encoded protein:

Check Swissmodel for a pre-existing close threaded model, else check out I-Tasser (best), Phyre (fast).

If you want a model of the DNA/RNA:

You are running blast without special settings —its not a short motif—, so these things may be long, so do be weary of any result. There are a few RNA folding prediction tools online (e.g. RNAFold).


Infer missing UTR features in GFF3 file I am working on a GFF file that is missing the 5'UTR and 3'UTR information. For example:

  ctg123 . gene      1050  9000  .  +  .  ID=gene00001;Name=EDEN
  ctg123 . mRNA      1050  9000  .  +  .  ID=mRNA00001;Parent=gene00001;Name=EDEN.1
  ctg123 . exon      1050  1500  .  +  .  ID=exon00002;Parent=mRNA00001
  ctg123 . exon      3000  3902  .  +  .  ID=exon00003;Parent=mRNA00001
  ctg123 . exon      5000  5500  .  +  .  ID=exon00004;Parent=mRNA00001
  ctg123 . exon      7000  9000  .  +  .  ID=exon00005;Parent=mRNA00001
  ctg123 . CDS       1201  1500  .  +  0  ID=cds00001;Parent=mRNA00001;Name=edenprotein.1
  ctg123 . CDS       3000  3902  .  +  0  ID=cds00001;Parent=mRNA00001;Name=edenprotein.1
  ctg123 . CDS       5000  5300  .  +  0  ID=cds00001;Parent=mRNA00001;Name=edenprotein.1

Is there a way to add rows in the GFF with the 5'UTR and 3'UTR ranges?

Among other things, the canon-gff3 program in the AEGeAn Toolkit infers and prints missing features like UTRs.

$ canon-gff3 < genes.gff3 > genes-with-utrs.gff3


Multiple correlation with R I have a data frame presenting different variables in column and for different species name in rows as follow. In the column a have only one independent variables(Proteome size) and all the others are independent variables. I want to use the PIC package to correlate all these dependent variables separately to my independent variables, a way to generate the R coefficient and p-value.

data <- read.table("t1.txt", row.names = 1, header = T,sep='')

data
path1 path2 path3 path4 Proteome_size
ahli 4.0 -3.2 0.4 -1.0 -0.2
alayoni 3.8 3.4 -1.8 2.2 -0.3
alfaroi 3.5 2.6 0.8 -2.4 -0.5
aliniger 4.0 0.1 -1.7 2.4 0.0
allisoni 4.4 2.0 -3.7 0.5 -0.1
allogus 4.0 -2.8 0.6 -1.0 0.6
altitudinalis 3.8 2.9 -6.1 2.3 1.2
alumina 3.6 0.7 1.5 -2.7 -1.3
alutaceus 3.6 1.2 -0.6 -1.7 0.3
angusticeps 3.8 4.6 -2.0 1.2 1.7

I have this code following code that can correlate any one of the ‘path*’ column to proteome size:

data <- read.csv("data.csv", row.names = 1
tree <- read.tree("anolis.phy")
Proteome <- data[, "Proteome_size"]
awe <- data[, "path1"]
names(Proteome_size) <- names(var1) <- rownames(data)
hPic <- pic(proteome_size, tree)
aPic <- pic(gene1, tree)
picModel <- lm(hPic ~ aPic - 1)
summary(picModel)

My question is that, how can I implement this code to do that for all the columns simultaneously and store the statistics info as (R and p_value) for each variable somewhere ?

Thanks

One option is to use one of the apply family functions. If you enclose your code above (which generates errors by the way) in a function, you can then "apply" this function on the rows or columns of your data with something like apply(X, MARGIN, FUN, …), the margin being 1 for rows or 2 for columns (as in your case).


What RNA-Seq expression value would be closest to Microarray equivalent? I know this question may seem strange.

I'm using Spearman correlation between gene expression profiles for various reasons (I won't go into details here). As a result, I often compare RNA-Seq and Microarray samples.
For preliminary analysis, I usually grab what version of the data is easily accessible (RPKM, FPKM,...), but I'd like to dig a bit more.

Intuitively, I'd think a value similar to RPKM or such would make more sense than raw counts, which is why I usually convert raw counts to RPKM (I know they're considered obsolete for statistical analysis of the data, what I'm doing relies on the correlation)

So I ask what RNA-Seq gene expression value would be closest to a microarray equivalent ?

i.e. What would (theoretically) maximize correlation between gene expression of the  same sample profiled with RNA-seq and microarray. RPKM ? FPKM ? TPM ? etc.

I think it is very hard to say which are the closest because they are not really comparable. But since you are using Spearman correlation, I guess RPKM, FPKM, and TPM do not change the order of gene expression levels. You might also want to normalize RNA-seq and microarray data so that they are more comparable.


getting KeyError on certain residues in BioPython: ATOM vs. HETATM lines When I want to load atoms for residuum (' ', 85, ' ') in chain A for protein 3tmm with structure[0][chain_id][residue], but I'm getting KeyError exception:

Traceback (most recent call last):
  [...]
  File "/home/pok/Checkouts/bilkoviny/loader.py", line 118, in _get_residue_atoms
    return list(self.structure[0][chain_id][residue])
  File "/usr/lib64/python2.7/site-packages/Bio/PDB/Chain.py", line 74, in __getitem__
    return Entity.__getitem__(self, id)
  File "/usr/lib64/python2.7/site-packages/Bio/PDB/Entity.py", line 40, in __getitem__
    return self.child_dict[id]
KeyError: (' ', 85, ' ')

I assume this is because this residuum do not have "ATOM" lines in PDB file:

[...]
ATOM    704  CD1 ILE A  84      19.348  -4.839  18.138  1.00 47.08           C
ANISOU  704  CD1 ILE A  84     5862   6108   5917   -216  -1513    609       C
HETATM  705  N   MSE A  85      21.619  -0.130  16.236  1.00 31.89           N
ANISOU  705  N   MSE A  85     5143   3072   3900   -766  -1011    962       N
HETATM  706  CA  MSE A  85      22.277   0.675  15.226  1.00 40.76           C
ANISOU  706  CA  MSE A  85     6708   3869   4911   -972   -729   1162       C
HETATM  707  C   MSE A  85      23.146   1.786  15.818  1.00 44.98           C
ANISOU  707  C   MSE A  85     7191   4141   5757  -1268   -558   1109       C
HETATM  708  O   MSE A  85      24.257   2.045  15.338  1.00 43.69           O
ANISOU  708  O   MSE A  85     7055   3855   5692  -1611   -151   1153       O
HETATM  709  CB  MSE A  85      21.244   1.265  14.269  1.00 46.50           C
ANISOU  709  CB  MSE A  85     8043   4309   5316   -709   -948   1345       C
HETATM  710  CG  MSE A  85      21.849   1.954  13.052  1.00 54.03           C
ANISOU  710  CG  MSE A  85     9604   4929   5995   -907   -623   1638       C
HETATM  711 SE   MSE A  85      20.503   2.755  11.877  1.00 84.35          SE
ANISOU  711 SE   MSE A  85    14348   8336   9366   -486  -1038   1884      SE
HETATM  712  CE  MSE A  85      21.690   3.538  10.518  1.00 83.79           C
ANISOU  712  CE  MSE A  85    15063   7865   8907   -915   -414   2326       C
ATOM    713  N   SER A  86      22.636   2.453  16.846  1.00 43.76           N
ANISOU  713  N   SER A  86     5234   6599   4794  -1778   -946    195       N
[...]

I have been told "MSE is same as MET, it is just with selenium instead of sulphur" and that I probably need to load these atoms (I'm approximating secondary structure element position, so I want individual atom coordinates).

Is there some way how to get atoms for this residue in this protein via BioPython?

You can find the answer in The Biopython Structural Bioinformatics FAQ:

What is a residue id?
This is a bit more complicated, due to the clumsy PDB format. A
residue id is a tuple with three elements:

The hetero-flag: this is 'H_' plus the name of the hetero-residue (e.g. 'H_GLC' in the case of a glucose molecule), or 'W' in the case
of a water molecule.
The sequence identifier in the chain, e.g. 100
The insertion code, e.g. 'A'. The insertion code is sometimes
used to preserve a certain desirable residue numbering scheme. A Ser 80
insertion mutant (inserted e.g. between a Thr 80 and an Asn 81
residue) could e.g. have sequence identifiers and insertion codes as
follows: Thr 80 A, Ser 80 B, Asn 81. In this way the residue numbering
scheme stays in tune with that of the wild type structure.

The id of the above glucose residue would thus be ('H_GLC', 100, 'A').

So in your case it would be ('H_MSE', 85, ' ').


What's the best way to download data from the SRA? Is it really this slow? I'm trying to download three WGS datasets from the SRA that are each between 60 and 100GB in size.  So far I've tried:

Fetching the .sra files directly from NCBI's ftp site
Fetching the .sra files directly using the aspera command line (ascp)
Using the SRA toolkit's fastqdump and samdump tools

It's excruciatingly slow.  I've had three fastqdump processes running in parallel now for approximately 18 hours.  They're running on a large AWS instance in the US east (Virginia) region, which I figure is about as close to NCBI as I can get.  In 18 hours they've downloaded a total of 33GB of data.  By my calculation that's ~500kb/s.  They do appear to still be running - the fastq files continue to grow and their timestamps continue to update.

At this rate it's going to take me days or weeks just to download the datasets.  Surely the SRA must be capable of moving data at higher rates that this?  I've also looked, and unfortunately the datasets I'm interested have not been mirrored out to ENA or the Japanese archive, so it looks like I'm stuck working with the SRA.

Is there a better way to fetch this data that wouldn't take multiple days?

Proximity to NCBI may not necessarily give you the fastest transfer speed. AWS may be deliberately throttling the Internet connection to limit the likelihood that people will use it for undesirable things. There's a chance that a home network might be faster, but you're likely to get the fastest connection to NCBI by using an academic system that is linked to NCBI via a research network.

Another possibility is using Aspera for downloads. This is unlikely to help if bandwidth is being throttled, but it might help if there's a bit of congestion through the regular methods:

https://www.ncbi.nlm.nih.gov/public/

NCBI also has an online book about best practises for downloading data from their servers.

On a related note, just in case someone sees this and EBI/ENA is an option, there's a great guide for how to do file transfer using Aspera on the EBI web site:

https://www.ebi.ac.uk/ena/browse/read-download#downloading_files_aspera

  Your command should look similar to this on Unix:

ascp -QT -l 300m -i <aspera connect installation directory>/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:<file or files to download> <download location>

In my case, I've just started downloading some files from a MinION sequencing run. The estimated completion time via standard FTP was 12 hours for about 32GB of data; ascp has reduced that estimated download time to about an hour. Here's the command I used for downloading:

ascp -QT -l 300m -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:/vol1/ERA932/ERA932268/oxfordnanopore_native/20160804_Mock.tar.gz .


After artificially creating events in a FASTA file, how do I keep track of the old coordinates? I'm beginning with the reference genome in FASTA format, hg19. I am reading the sequence into a Python dictionary with BioPython:

genome_dictionary = {}
for seq_record.id in SeqIO.parse(input_fasta_file, "fasta"):
    genome_dictionary[seq_record.id] = seq_record

This creates a Python dictionary with keys as chromosome names and values as strings from the FASTA.

My goal is to take the FASTA, and simulate artificial SNPs and InDels along the FASTA to benchmark various algorithms. If I was only creating SNPs in the FASTA, this would be a trivial problem. The reference FASTA has the original coordinates, and the SNPs would also be in a coordinate system which is the same.

However, when I randomly insert an insertion (e.g. an insertion of 12 bp), I move all events into a new coordinate system such that

new coordinate system = old coordinate system + 15

This becomes very difficult to keep track off throughout the entire genome after several interations. I therefore cannot keep track of the new and old coordinate system without running an algorithm and creating a VCF, which may have errors. This defeats the purpose of the simulation.

The output would be either a BED file or VCF which has keep track of the changes.

What method/data structure is used to solve this problem? Would LiftOver work? Maybe CrossMap?

http://genome.sph.umich.edu/wiki/LiftOver
http://crossmap.sourceforge.net/

I've written a handful of programs from scratch to simulate mutations and variations in real or simulated sequences.

The trick has always been to sort the variants by genomic coordinate, apply the variant with the largest coordinate first, then apply the variant with the second largest coordinate, all the way down to the variant with the smallest coordinate.

Indels and other variants that affect sequence structure only affect subsequent coordinates. So going in reverse order ensures variants at the beginning of the sequence don't throw off variants at the end of the sequence.


Does Prokka do six-frame translations? Does Prokka do six-frame genome annotations? The corresponding publication does not mention it.

  Seemann, Torsten. 2014. “Prokka: Rapid Prokaryotic Genome Annotation.” Bioinformatics  30 (14): 2068–69.

Yes, Prokka will consider all six-frames, the gene prediction is actually done with Prodigal, paper:
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2848648/

github:
https://github.com/hyattpd/Prodigal


High Phred Quality score VCF file I am working on ION Torrent S5 XL and I am facing to an issue with VCF output generated by the variant caller of thermofisher. It stating on VCF 4.1 that Phred Quality Score is on a range of 0 to 255. In my VCF I have values over than 255. Do you know if there is a way to transform these values?

See this post on the GATK forums.

It is normal to see values above 255.  The QUAL score is the posterior probability that all samples in your callset are homozygous reference. I quote Sheila's response from the thread posted above:

  The QUAL score is the Phred-scaled posterior of AC = 0. We use the AC priors and the PLs to get the likelihood of the data given each AC, then use those to get the posterior probability for each AC. From there, the calculation is 1 - Pr{AC > 0}.

Please see the GATK post in Methods and Algorithms for the nitty gritty details. Beware, Bayesian statistics abound.


Sequence allignment with suffix array? I'm writing the code from scratch and not using libraries for the actual indexing/search since it's a project for school. Any blackbox advice is fairly useless.
I'm planning on aligning some short reads with the D. melanogaster genome but, since I was thrown into this project with very little info/ experience I'm stuck on a few things.
Firstly, for preprocessing the fasta file:

What should I do with Ns (ambiguous bases), I've seen some people removing them or replacing them by gaps but this does not seem like the same information? Should they be left in?

As far as I understand, the lowercase letters stand for potential repeats. However, If I just ignore this, it is possible my sequences might not align to the repeat because it is not actually stated as a repeat. E.g: "ATATATAT" will not align to 'at'.

As far as I understand the difference between using suffix array and suffix tree: suffix tree will search in linear time (but is memory intensive), while suffix array will be space efficient (but searches through ALL sequences each time). I'm leaning towards suffix array because I'm working on my laptop and memory is kind of limited. I'm also thinking that, since you sort the array lexicographically, you really only have to check the sequences if they start with the same first base as the read you're searching for?
Another question on the arrays: if my genome has 100M base pairs, the array will have 100M entries, some of them being close to 100M bases long? Should I be using the seed and extend method?
On the other hand, since the fruit fly genome is one of the smaller genomes, maybe it would be feasible to use a tree?
Is there any way to validate the alignment quickly at the end? How much error should I be allowing?
This were a lot of question, I've been trying to read up on things but I can mostly only find applications using biopython etc.
Any help would be hugely appreciated!

Suffix arrays vs Suffix tree is a problem quite often discussed in bioinformatics classes, but it's not that much used in practice. Most of people use already optimized mappers bowtie2 or bwa-mem.
I would advise you against your intention of mapping the reads using your own implementation of anything. You don't need to reinvent the wheel, just use an existing solution.
By the way, the Ns in the reference are either hard-masked nucleotides (that would be repeats of some sorts) or (more often) simply unknown nucleotides. So there is no way to map your reads there anyway (and the mentioned mappers will handle the Ns for you).


Highly heterozygous reads mapping I have short (67bp) Hi-C reads from a highly heterozygous organism (~15% between-haplotype SNP divergence) and I have both reference haplotypes.

I wanted to try and compare different haplotyping softwares for Hi-C reads using these reads as benchmarking dataset. When mapping the reads separately to each haplotype, I get good mapping statistics. When I map the reads to a single haplotype with all the heterozygous SNPs masked (into N's), I get  very poor mapping rates.

I would like to be able to map the reads when the real haplotypes are unknown (reference is a mix of haplotypes).

I use minimap2 to map the reads with the sr preset. I tried decreasing the mismatch penalty (-B) to 1 and increasing the --score-N value, but this had no effect.

As shown in the attached IGV screenshot, coverage drops to 0 when the SNP density increases.
Is it feasible to map reads with such a high heterozygosity on a single (masked) reference ? Should I use another tool ?

I believe you may be able to map your reads, but I don't know how to do that with minimap2.

I recommend running gsnap, which is more SNP tolerant and provides a number of parameters which are likely to help.

For instance, I believe that most aligners will treat 'N' characters as mismatches when you align. GSNAP has parameters to account for this.

--query-unk-mismatch=INT       Whether to count unknown (N) characters in the query as a mismatch
                               (0=no (default), 1=yes)
--genome-unk-mismatch=INT      Whether to count unknown (N) characters in the genome as a mismatch
                               (0=no, 1=yes (default))

It also has a mismatch parameter similar to the one you described for minimap2.

-m, --max-mismatches=FLOAT     Maximum number of mismatches allowed

Try running with the 'genome-unk-mismatch' parameter above (with the masked reference). That may be your best bet. There are also other parameters that may help, but this should be a good start.


How can I calculate loss of heterozygosity (LOH) in NGS sequencing data? I'm analyzing a tumor sample and a healthy sample from the same patient. I want to use sciClone to look at tumor clonality. One input is the genomic regions that need to be excluded due to LOH. I have both WES and WGS.

Is there currently a "gold standard" for which tool will return the most robust regions?

The sciClone paper mentions VarScan2, APOLLOH, HMMScan? Are these still state of the art?

Is it feasible to calculate LOH programmatically without the need of external software?

Answer from @fabio-marroni, converted from comment:
WARNING! Self-promotion! We published a freely-available tool for detection of LOH (and also redeuction fo homozygosity). Here you find the link to the paper and therein will find the link to the software. We compared with other tools that you might also try to see which one is the best performing. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5466433/


Apply trained PRS on another dataset I am using PRSice to compute the PRS over a train set and want to use the coefficient used on the train set to apply it on another set which I will call the test set.
Once I compute the PRS I get a file called analysis.prsice
This is the few first lines

Pheno Set Threshold R2 P Coefficient Standard.Error Num_SNP

- Base 5e-08 0.00550157 0.0659631 -21.0513 11.4492 43

- Base 5.005e-05 0.00396637 0.113759 -61.8545 39.1108 184

- Base 0.00010005 0.00236122 0.220463 -68.9527 56.2742 286

- Base 0.00015005 0.00113035 0.394503 -59.4058 69.7677 382

- Base 0.00020005 0.00072253 0.495129 -55.2605 81.0064 471

- Base 0.00025005 0.000334179 0.642121 -41.5674 89.4432 538

- Base 0.00030005 7.5306e-05 0.825111 -21.391 96.8021 601

- Base 0.00035005 3.95773e-06 0.959527 5.25804 103.612 668

- Base 0.00040005 0.000175506 0.735018 37.723 111.455 743

I thought that I could apply the PRS on another dataset by multiplying the coefficient in the result by the allelic dosage of the corresponding SNPs. However, I am not sure how to know which are the corresponding SNPs of the coefficients.
Also, if there was a command that does that and it would be easier, could you please tell me how to do it?
Thanks

The coefficient is the third column from the right, you can tell because some are positive and other negative. When the value is negative the association is a negative correlation. When it is positive (last two) that is more like the result you are seeking. The further away from zero the stronger the correlation (positive or negative)
These are very large values for the coefficients. The SNP is the last column on the right. I would disgard the top two results because the probabilities are too borderline, e.g. if you used Bonferoni's correction.
There isn't enough information on this calculation nor the data set to make any further insights other than you have 7 SNPs showing PRS


Comparing AUCs: Discrimination of same Control from different Test Group - paired or unpaired? miRNA I want to compare two AUCs using a bootstraping method from pROC package (roc.test).
I have done a classification using Logistic Regression and Serum miRNAs as prediction factors. The Controls are in both cases the same, the test group is different.
1. Cancer (high risk) vs. Control = AUC 1
2. Cancer (low risk) vs. Control = AUC 2

If I want to compare this two AUCs - should I use a paired or unpaired test?
I would tend to using a paired test, because of exact the same Control group (Control groups are correlated).

Generally its paired. The bootstrapping personally I find concerning. I presume you are using it to augment your data set.
Bootstrapping at best tends to be a controversial statistical operation and could result in the non-significance you are observing. Anyway if its not significant between high and low - thats your answer.
You could of course perform a single ML classifying low, and high cancer risk against control. That would provide a method to understand the associations of miRNA.

Bootstrapping issue, I had a think about this. It depends on how the bootstrapping is done. In a paired sample this (I assume) samples with replacement at each paired position for both ROCs. If this is not the case - the test is junk. However, if this is the case, its actually a good test, because the paired structure of the data is always preserved no matter how many bootstrap samples are performed. So like if a 100 were done and the pairwise statistic (which we don't know what it is - but something like KS test) and its >0.05 for all 100 samples - I don't think there's any doubt thats a robust result.
Again if the two ROCs are bootstrapped independently, in my opinion the test is junk.
Reason The concern of bootstrapping is it causes weird distributions which were not present in the original data. This is only okay if there's lots of experience about interpreting the result. For example, in phylogenetics there's decades of experience interpreting bootstrapping and a consensus about what it means. This often isn't the case for most other areas of statistics.
I note this particular issue is really ramping up the rep score ;0)