How do I set a neural network to loop multiple times and average the resulting values? I have a script in R/RStudio which creates random datasets of binomial variables, feeds them through a neural network, and calculates their likelihood ratio statistic and deviance. I'd like the script to loop with the seed incremented by 1, and collect the value from each run into a data frame, from which the values get averaged. For example, if the script returns 5, 10, and 15 on the first, second, and third iterations of the loop, I'd like the script to create a data frame with those values and then calculate their mean.
Here is my script, annotated:
library(neuralnet)

x=1021

#change nrow to number of times code is looped
result <- data.frame(matrix(nrow = 1000, ncol = 2))
colnames(result) <- c("LogL", "D")

repeat{

set.seed(x)

#2 inputs, max number of terms 3
x1=rbinom(2000,1,0.5)
x2=rbinom(2000,1,0.5)
y=rbinom(2000,1,0.5)

data=data.frame(x1,x2,y)

head(data)

#hidden units = 2
nnet=neuralnet(y~x1+x2,data,hidden=2,err.fct = "ce",linear.output = F,likelihood = T)

predictions=compute(nnet,data[,1:2])

predictions$net.result

#likelihood ratio stat
argument=y*log(predictions$net.result)+(1-y)*log(1-predictions$net.result)

head(argument)

#sum of likelihood ratio stat
LogL=sum(argument)

LogL

D=-2*LogL

D

result[x-1022, 1] <- LogL
result[x-1022, 2] <- D

result

#results <- data.frame(LogL, D)

x=x+1
if (x==2021) {
  break
}

print(x)

}

result

My repeat loop works, but a thousand iterations, for example, takes multiple minutes. Is there a more efficient way to code such a loop?

The concern with your approach is called 'leakage' because you are parameterising the same data set that you are training. This can easily lead to overtraining is in ANN (artificial neural networks) and is not cool. Overtraining is where the ANN very tightly wraps around the training data set that it performs poorly on any other data set. In other words it is overtightened.
The way they get around it is to split the data set usually 60:20:20 into a training split, parameterisation (validation) split and test (testing) split. If the parameterisation split gets mixed with the training set - leakage happens and overtraining results. Normally the data is trained, and the trained model parameters are assessed/optimised in the parameter split. The trained data and optimised parameters are then tested in the test split.
As long as you are aware of the approach you are using, that I guess that is ok. If you separate your training and parameterisation (validation) data that would be cool.


DMPfold: is it suitable to get 5000 models to compare with rosetta I ran this command for 10 iteraions and 500 models.

run_dmpfold.sh sequence.fasta sequence.21c sequence.map ./dmpfold-lysin 10 500 &> dmpfold-log.log

I am looking to generate 5000 dmpfold models to compare rosetta results with dmpfold; is it reasonable to do this? and, 5000 models can be generated by using "10 5000" parameter or "10 500" when running dmpfold as stated in the command line? Thanks!

DMPfold tends to produce low-variability models - it produces a small number of distinct conformations, and either the right one is in there or it isn't. This arises from the way the models are generated from distance constraints. Rosetta, by contrast, produces many different conformations which is why 5000 or more models tend to be generated. This arises from Rosetta's method of assembling small fragments, which allows for many more combinations.

You could generate 5000 models with DMPfold. Running with 10 5000 would give you 5000 models as output from the last iteration. Running with 10 500 would give you 500 models from the last iteration but you could also use the 500 models from each of the previous iterations (saved as ensemble.1.pdb etc.) to get 5000 models overall.

However I wouldn't recommend it, it sounds like a waste of disk space and compute time. Just run with default parameters 3 50 and use the final models saved as final_*.pdb.

For some numbers on how DMPfold compares to larger numbers of Rosetta models, and how much better DMPfold gets when you use more models, see Table 1 of the paper.


SNPEFF genome build with specific codon usage Failing to build my own Arabidopsis thaliana reference with Mt using 'codon.Mitochondrial' and Pt using 'codon.Bacterial_and_Plant_Plastid' but the same commands works well with both lines with codon tables removed below.

The snpeff commands

java -Xms24g -Xmx24g -jar /opt/biotools/snpEff/snpEff.jar build -gff3 -v Arabidopsis_thaliana.TAIR10.46

returns:

java.lang.RuntimeException: Error parsing property 'Arabidopsis_thaliana.TAIR10.46.Pt.codonTable'. No such codon table 'codon.Bacterial_and_Plant_Plastid'
        at org.snpeff.snpEffect.Config.createCodonTables(Config.java:169)
        at org.snpeff.snpEffect.Config.readConfig(Config.java:650)
        at org.snpeff.snpEffect.Config.init(Config.java:480)
        at org.snpeff.snpEffect.Config.<init>(Config.java:117)
        at org.snpeff.SnpEff.loadConfig(SnpEff.java:451)
        at org.snpeff.snpEffect.commandLine.SnpEffCmdEff.run(SnpEffCmdEff.java:1000)
        at org.snpeff.snpEffect.commandLine.SnpEffCmdEff.run(SnpEffCmdEff.java:984)
        at org.snpeff.SnpEff.run(SnpEff.java:1183)
        at org.snpeff.SnpEff.main(SnpEff.java:162)

(same error for the Mt when removing Pt)

my fasta has the two chromosomes

>1 dna:chromosome chromosome:TAIR10:1:1:30427671:1 REF
>2 dna:chromosome chromosome:TAIR10:2:1:19698289:1 REF
>3 dna:chromosome chromosome:TAIR10:3:1:23459830:1 REF
>4 dna:chromosome chromosome:TAIR10:4:1:18585056:1 REF
>5 dna:chromosome chromosome:TAIR10:5:1:26975502:1 REF
>Mt dna:chromosome chromosome:TAIR10:Mt:1:366924:1 REF
>Pt dna:chromosome chromosome:TAIR10:Pt:1:154478:1 REF

my data folder is OK

data/Arabidopsis_thaliana.TAIR10.46/
drwxr-xr-x   3 u0002316 domain users 4.0K Apr 15 13:33 .
drwxr-xr-x 244 u0002316 domain users  16K Apr 10 09:39 ..
-rw-r--r--   1 u0002316 domain users 106M Apr 10 09:40 genes.gff
-rw-r--r--   1 u0002316 domain users 117M Apr 10 09:40 sequences.fa

my snpeff.config contains:

Arabidopsis_thaliana.TAIR10.46.genome : Arabidopsis_thaliana.TAIR10.46
Arabidopsis_thaliana.TAIR10.46.Mt.codonTable : codon.Mitochondrial
Arabidopsis_thaliana.TAIR10.46.Pt.codonTable : codon.Bacterial_and_Plant_Plastid

the two codon tables are present by default in the snpeff.config

codon.Mitochondrial                         : TTT/F, TTC/F, TTA/L+, TTG/L+, TCT/S, TCC/S, TCA/S, TCG/S, TAT/Y, TAC/Y, TAA/*, TAG/*, TGT/C, TGC/C, TGA/W, TGG/W, CTT/L, CTC/L, CTA/L, CTG/L+, CCT/P, CCC/P, CCA/P, CCG/P, CAT/H, CAC/H, CAA/Q, CAG/Q, CGT/R, CGC/R, CGA/R, CGG/R, ATT/I+, ATC/I+, ATA/I+, ATG/M+, ACT/T, ACC/T, ACA/T, ACG/T, AAT/N, AAC/N, AAA/K, AAG/K, AGT/S, AGC/S, AGA/R, AGG/R, GTT/V, GTC/V, GTA/V, GTG/V+, GCT/A, GCC/A, GCA/A, GCG/A, GAT/D, GAC/D, GAA/E, GAG/E, GGT/G, GGC/G, GGA/G, GGG/G

codon.Bacterial_and_Plant_Plastid           : TTT/F, TTC/F, TTA/L, TTG/L+, TCT/S, TCC/S, TCA/S, TCG/S, TAT/Y, TAC/Y, TAA/*, TAG/*, TGT/C, TGC/C, TGA/*, TGG/W, CTT/L, CTC/L, CTA/L, CTG/L+, CCT/P, CCC/P, CCA/P, CCG/P, CAT/H, CAC/H, CAA/Q, CAG/Q, CGT/R, CGC/R, CGA/R, CGG/R, ATT/I+, ATC/I+, ATA/I+, ATG/M+, ACT/T, ACC/T, ACA/T, ACG/T, AAT/N, AAC/N, AAA/K, AAG/K, AGT/S, AGC/S, AGA/R, AGG/R, GTT/V, GTC/V, GTA/V, GTG/V+, GCT/A, GCC/A, GCA/A, GCG/A, GAT/D, GAC/D, GAA/E, GAG/E, GGT/G, GGC/G, GGA/G, GGG/G

Anybody ever succeeded to build a custom genome with specific codon tables?

Thanks in advance

codon table issue

I solved the issue with codon tables, the correct way to specify a table is the following:

Arabidopsis_thaliana.TAIR10.46.genome : Arabidopsis_thaliana.TAIR10.46
Arabidopsis_thaliana.TAIR10.46.Mt.codonTable : Mitochondrial
Arabidopsis_thaliana.TAIR10.46.Pt.codonTable : Bacterial_and_Plant_Plastid

Note that one should remove the leading prefix 'codon.' from the name of the tables chosen from the config file

Concerning accessory files:

data to reproduce this issue:

baseurl="ftp://ftp.ensemblgenomes.org/pub/plants/release-46"
wget ${baseurl}/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
wget ${baseurl}/fasta/arabidopsis_thaliana/cds/Arabidopsis_thaliana.TAIR10.cds.all.fa.gz
wget ${baseurl}/fasta/arabidopsis_thaliana/pep/Arabidopsis_thaliana.TAIR10.pep.all.fa.gz
wget ${baseurl}/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.46.gff3.gz
wget ${baseurl}/variation/vcf/arabidopsis_thaliana/arabidopsis_thaliana_incl_consequences.vcf.gz

I suspect that the format used by Ensembl is not 100% OK for SnpEff; I therefore simplified the fasta headers to keep only the first argument (name) and remove the remaining of the line

eg.

>AT3G05780.1 pep chromosome:TAIR10:3:1714941:1719608:-1 gene:AT3G05780 transcript:AT3G05780.1 gene_biotype:protein_coding transcript_biotype:protein_coding gene_symbol:LON3 description:Lon protease homolog 3, mitochondrial [Source:UniProtKB/Swiss-Prot;Acc:Q9M9L8]

becomes

>AT3G05780.1

Despite this modification, I still have an issue with

java -jar /opt/biotools/snpEff/SnpEff.jar build -gff3 -v Arabidopsis_thaliana_TAIR10_46 2>&1 | tee -a log.txt

a LOT of warning about not being able to add UTR while it seems that UTR's are well present in the GFF3 file

  Reading GFF3 data file  :
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  WARNING: Cannot find transcript 'transcript:AT1G01010.1'. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  18    '1  araport11   five_prime_UTR  3631    3759    .   +   .   Parent=transcript:AT1G01010.1'
  WARNING: Could not add UTR. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  18    '1  araport11   five_prime_UTR  3631    3759    .   +   .   Parent=transcript:AT1G01010.1'
  WARNING: Cannot find transcript 'transcript:AT1G01020.3'. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  69    '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.3'
  WARNING: Could not add UTR. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  69    '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.3'
  WARNING: Cannot find transcript 'transcript:AT1G01020.5'. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  88    '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.5'
  WARNING: Could not add UTR. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  88    '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.5'
  WARNING: Cannot find transcript 'transcript:AT1G01020.4'. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  109   '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.4'
  WARNING: Could not add UTR. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  109   '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.4'
  WARNING: Cannot find transcript 'transcript:AT1G01020.1'. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  128   '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.1'
  WARNING: Could not add UTR. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line
  128   '1  araport11   three_prime_UTR 6788    6914    .   -   .   Parent=transcript:AT1G01020.1'
  WARNING: Cannot find transcript 'transcript:AT1G01030.1'. File
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/genes.gff'
  line 155  '1  araport11   three_prime_UTR 11649   11863   .   -   .   Parent

nothing good coming out of the cds.fa input although I treated it with the same name fix as the pep file

  00:00:19  CDS check (FASTA file):
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/cds.fa'

  00:00:19  Checking database using CDS sequences 00:00:19  Reading CDSs
  from file
  '/opt/biotools/snpEff_v4.4/./data/Arabidopsis_thaliana_TAIR10_46/cds.fa'...
  00:00:20  done (75976 CDSs). 00:00:20 Comparing CDS...    Labels:         '+' :
  OK        '.' : Missing       '*' : Error
    ....................................................................................................
    ....................................................................................................
    ....................................................................................................

  many more lines with dots
    ....................................................................................................
    ....................................................................................................
    ....................................................................................................
    ...............................     CDS
  check:    Arabidopsis_thaliana_TAIR10_46  OK: 0   Warnings: 0 Not found:
  104147    Errors: 0   Error percentage: NaN% Fatal error: No CDS checked.
  This is might be caused by differences in FASTA file transcript IDs
  respect to database's transcript's IDs. Transcript IDs from database
  (sample):     'transcript:AT3G12250.7'    'transcript:AT5G18570.1'
    'AT3G27810.1'   'AT3G27810.2'   'transcript:AT2G16810.1'
    'transcript:AT3G01170.1'    'AT4G37030.1'   'transcript:AT1G79940.4'
    'transcript:AT2G44850.1'    'AT1G49540.1'   'AT1G49540.2'
    'transcript:AT4G35290.5'    'transcript:AT4G15840.2'
    'transcript:AT4G27500.1'    'AT4G37022.1'   'transcript:AT2G40008.1'
    'AT3G27820.1'   'transcript:AT4G24265.1'    'transcript:AT2G20618.1'
    'transcript:AT4G18375.2'    'AT4G37020.2'   'AT4G37020.1' Transcript IDs
  from database (fasta file):   'AT4G02541'     'AT1G12410.1'   'AT4G26510'
    'AT4G02540'     'AT3G04500.1'   'AT3G04500.2'   'AT3G04500.3'   'AT2G36250'
    'AT5G63350'     'AT2G36255'     'AT2G12280'     'AT1G03970.1'   'AT5G18390.1'
    'AT3G04490.4'   'AT5G38386'     'AT3G04490.2'   'AT3G04490.3'
    'AT3G04490.1'   'AT5G38380'     'AT4G02530'     'AT4G26520'     'AT1G69280.2'

link to the full stdout log after running the command without cds.fa

https://www.dropbox.com/s/gxa4e533p64wxm9/SnpEff_build-log.txt?dl=1


In Uniprot how do I know which variants are mutations and which are polymorphic? In Uniprot how do I know which variants are mutations and which are polymorphic?
I am looking at P38398 (BRCA1_HUMAN) specifically.

Many thanks.

This answer is very late but may still help others:
Have a look at the UniProt index of human polymorphisms and disease mutations at https://www.uniprot.org/docs/humsavar

Please don't hesitate to contact the UniProt helpdesk if you have any additional questions.


How to deal with zeros in the denominator when calculating percent enrichment? I have a data.frame with results from a CRISPR screen.

I want to calculate the % enrichment of a variable of interest var from baseline to treatment group.

Here is the distribution of var in linear:

I can do scale var with log transformation

But many rows with var = 0 will be lost.
Moreover, I cannot compare % enrichment by dividing treatmentgroup::var to Control/baseline::var because I will run into the problem of dividing by zero.

I've seen people circumvent this problem by adding a small value ex)0.01 to the entire column of var.
But I'm wondering if there's a better way of circumventing this problem. For example, a simple transformation that can deal with zeros and negative values.

Appreciate any input!

Cube root transformation is another way to deal with zeros without dropping them.

var^(1/3)

EDIT as Michael pointed out, square root transformation is another--perhaps more common--transformation that leaves zeros as zeros.

sqrt(var)


Introduce errors in reference transcripts according to external dataset error model I would like to modify some reference transcripts from Ensembl (D. melanogaster) to introduce a controlled rate of random errors in the sequences. The idea would be to introduce random base substitutions in these sequences, no indels for now, because I would like to keep the transcript sequence length as it is in the reference.

The rate of error per transcript will be determined according to an error profile computed from an external set of RNA-seq reads (e.g., generated with ONT MinION)

The aim of this modification would be to establish a rough benchmark of the performances of aligners to use over transcripts from spliced reads (rna-to-genome), aka with more than one exon.

Any idea of which software would be best for this purpose?

Do any of the answers for this question help? Karel Brinda has mentioned a few read simulators in the answer to that question, and has a thesis with more information.

Excluding INDEL errors doesn't sound like a good idea; length can still be preserved even if doing that, it just needs an adjustment at the end of the sequence. Note that if you're trying to model nanopore reads, what you're really modelling is the base-caller, rather than the sequencer. I mention this in more detail in my answer.

In most cases where errors are modelled, I find it better to use publicly-available data instead. Especially for nanopore data, there are unmodelled systematic errors in the base-callers and sequencer that can't be simulated using any programs (because they are unmodelled). The following paper would be a good place to start for cDNA sequences, which looks at single-cell data from mouse (C57Bl/6) B1a cells:

http://www.biorxiv.org/content/early/2017/04/13/126847

Illumina and ONT reads for that study can be found in SRA under accession number SRP082530.

I don't know of any recent D. melanogaster studies that have been done using nanopore. There's always the option of spending $1000 on a purchase of a MinION with an RNA starter kit to do the study yourself. Here's an older targeted gene study, but bear in mind that it was using an R7.3 flow cell, so errors rates will be much higher than what is currently available:

https://genomebiology.biomedcentral.com/articles/10.1186/s13059-015-0777-z


Should biological replicates be the most similar pairs in an RNAseq experiment? The attached figure (from deseq2) shows the sample to sample distances of a RNAseq experiment with 4 conditions (A,B,C,D) at 2 time points (0h,4h) with 2 biological replicates.

I am a little bit puzzled that the biological replicates are not always the most similar pairs. In a perfect world biological replicates should be identical but I assumed that at least they should be very similar.

Is this something I have to worry about, or is this more or less normal?

In summary you basically have noisy data and to simply answer your question it depends on the genetic diversity of the samples you are dealing with. Note, from my recollection this was confusing to start with because the "replicates" column was labeled as "treatment", but I could be wrong.

Revised answer, firstly there is zero stationarity - thus its a straight time-series.

Summary As a time series

The statistical noise between A and D is lower than the difference in phenotypic expression between each and all othe replicates in the analysis ABCD - FOR DAY 0 VERSUS DAY
4, so they cluster and are distinct.
The statistical noise between B and C is greater than the phenotypic variation between their expression profiles for DAY 0 to DAY 4. I suspect they are genetically very similar. Thats the conclusion, you can spice it up with different analysis, but I'd leave it because might find a lack of statistical significance in your core observations (A and D).
Within the time-points the data is very noisy and doesn't appear to show a distinct signal. Thats okay, this is essentially a time series.
Sample D is distinct in all scenarios (which is good)

Result To reiterate and expand a bit, the result you have is that sample A and D biological replicates are consisent within the analysis you are performing, i.e. day 0 verses day 4 (below). Thus they cluster together and this is proof of distinct phenotypes above the statistical noise (A with A and D with D). An intra-time point dendrogram would have a different result. The replicates samples B and C are not consistent and the interpretation is that the phenotypic expression profiles are so similar they are essentially identical within the boundaries of statistical noise.

Analysis The dendrogram is combines day 0 to day 4 analysis with day 0 to day 0 analysis and day 4 to day 4 analysis. HOWEVER, the clustering algorithm will correct for this and separate based on the largest distance (day 0 versus day 4). The signal that you get in your dendrogram is the largest signal in the data set which is day 0 vs day 4, and this is how UPGMA clustering works (I probably have previous posts on the statistical algorithm of UPGMA). So actually that is fine.

Advice With data this noisy I'd go no further in attempting to deploy critical probablity (a p value)  because this data could be the result of a complex distribution. This would require a statistician. Personally
 I'd upload the PCA as a separate post, but its just my opinion, because UPGMA versus eigenvalue based PCA are statistically quite different approaches.


bcftools output vs. bgzip I've got what seems like it should be a simple question, but I can't seem to figure it out from google. bcftools has 4 output options:

Output compressed BCF (b), uncompressed BCF (u), compressed VCF (z), uncompressed VCF (v).

Are any of these the same as the output from bgzip?

Yes. I haven't found it documented anywhere, but I ran a simple test:

$ bcftools view -Oz --no-version foo.vcf > bar.vcf.gz
$ md5sum bar.vcf.gz
18a2f494710cf170ba79892936015ba8  bar.vcf.gz

$ gunzip bar.vcf.gz
$ bgzip bar.vcf
$ md5sum bar.vcf.gz
18a2f494710cf170ba79892936015ba8  bar.vcf.gz

$ gunzip bar.vcf.gz
$ gzip bar.vcf
$ md5sum bar.vcf.gz
9f5b1cf654966912b78a9cb7175ffc2c  bar.vcf.gz

I passed a vcf file through bcftools view -Oz, checked the md5sum of the resulting file, then uncompressed it and recompressed with bgzip and confirmed that the resulting file had the same md5sum. Finally, I also compressed using normal gzip to be sure that the md5sum would be different, as it was.

So yes, it does indeed look like bcftools creates files that are the same as though you'd passed them through bgzip.


QC measures for NGS sequencing What are good means for performing quality control (QC) or NGS reads?

I'm aware of:

FastQC
NGS Screen
Kraken (e.g., for screening against contaminants)

What are other popular means for such QC?

MultiQC can merge all your different reports into a single one.
Which could be useful once you manage to know which QC tools to use.


Loading all PPI information from Uniprot into a graph in NetworkX I am looking for guidance to load all SwissProt PPI data for Homo sapiens into a graph (from NetworkX). In particular, the part that I'm struggling with is gathering the data from SwissProt for each human protein.

The data that I'm looking to extract is protein-protein interaction. This is not included in the files found here, which @rtviii (answer below) suggested to use. Rather, when say you click on an entry in SwissProt, for each entry, there is a section devoted to "Interaction" as shown in the image attached below. I would like to extract that information for each human protein, and load it into a graph data structure using NetworkX

EDIT 2:

How I ended up doing it:

    query = "reviewed:yes+AND+organism:9606"
    args = {'query': query, 'format': 'xls', 'columns':'id,entry name,interactor'}

    uniprot_url = "https://www.uniprot.org/uniprot/?{}".format(urllib.parse.urlencode(args))

    with open(<DESIRED FILE PATH>, "wb") as target:
          target.write(requests.get(url).content)

The query link for all interactors would be thus : https://www.uniprot.org/uniprot/?query=organism%3A9606+AND+interactor%3A*&sort=score

There is always an option of downloading the static version of the service if you find it incovnient to deal with their API for your use case from here and have it all to yourself, pipe data into your dataframes locally. I find it useful to just broswe their repositories sometimes instead of reading the specifications.

In general, you can get a variety of datashapes  from uniprot that you cannot get from their static database-image by specifying the appropriate parameters in the request. You can browse and explore parameters by expanding their GUI's search fields as described in this little guide. Then you can request your data like so, in python for example.

import requests
URL = "https://www.uniprot.org/uniprot/"
PARAMS = {'organism': 9606, 'format': 'tab', ....#Additional parameters#...}
response = requests.get(url=URL, params=PARAMS)


What is the etymology of "Entrez ID"? Since I have seen NCBI gene names were called "Entrez ID" for the first time, I am wondering where that comes from. Such a weird name!

Does anybody know where that originates?

My hypothesis is: in French, "entrez" can be translated to "please enter" (it is to enter in the imperative mood), so I'm guessing that might come from a past website/database in French where Entrez ID as in "Please enter the ID" was written in the search bar...

But I could not find any reference... I cannot be the only one to ask himself this!

Update

As pointed out by neilfws:

According to Wikipedia and Practical bioinformatics by Michael Agostino: The name "Entrez" means "Come in!" in French (pronounced “on-tray”); it was chosen to reflect the spirit of welcoming the public to search the content.

Old (probably incorrect) theory

The first version of Entrez database was distributed by NCBI in 1991 (on CD-ROM). At that time, it consisted of nucleotide sequences, protein sequences and associated citations and abstracts from MEDLINE.

So even at the early stage it contained objects (entries) of very different types (sequences and abstracts). It is hard to find a word which will accurately and briefly describe all the types of data.

While I don't know the correct answer and this is purely a speculation, my guess is that Entrez is an intentionally misspelled word "entries"; which accurately describes all the datatypes in the database.

Using a misspelled word as a name is not uncommon, e.g. Google comes from misspelled "googol" and Flickr is a misspelled version of "flicker".


Most scalable pseudotime ordering algorithm What algorithms for linear Pseudotime trajectory construction (diffusion-based) are the most scalable to large datasets?

I'm currently using Slingshot based on the recommendation in this manuscript: https://www.biorxiv.org/content/biorxiv/early/2018/03/05/276907.full.pdf

Slingshot is nice and easy to use, but it turns out it is terrible at scaling to big datasets.

For 2000 highest variable genes and a random sampling of a variable number of cells in my ~80% sparse expression matrix I get the following computation times in R (96 GB RAM):

#cells / Time (min)
100 / 0.7
1000 / 4.3
2000 / 12.3
5000 / 71.1

At this rate it would take nearly 1.5 days to process 200,000 cells and 4.5 months to process 2 million cells based on 2,000 HVG.

I feel like Slingshot is limited by a fundamentally unnecessary comparison of diffusion for all cells against all cells, rather than preclustering and then finer resolution of pseudotime across cluster edges.

What algorithms are better, but still implement a similar realization of pseudotime ordering?  I just need a single linear trajectory.

As it turns out, pseudotime algorithms are exceedingly difficult to parallelize. Pseudotime does not scale linearly to a subset of data compared to the full dataset, and as such a course-grained approach to parallelizing pseudotime calculations on random bins of the data, followed by bin merging, does not work. These values will never converge on the actual full-dataset solution.

I put in a pretty good effort trying to parallelize slingshot, but with no success. I realize my code may be of limited utility without more explanation, but I'm not putting more time into this. FWIW:

This code:

Finds computationally optimal bin size for slingshot based on a number of genes and features
Randomly splits the SingleCellExperiment object into that number of bins and runs slingshot on each bin
Orients each bin so that a cell with a pseudotime of 1 in the first bin is similar to a pseudotime of 1 in every bin. Pseudotime can get flipped depending on the sampled subset.
Merges all bins into a single matrix so each cell has a assigned pseudotime based on randomly sampled bins.
Repeats steps #2 to #5 with a different random seed for binning every time.
Measures convergence towards actual solution at three pseudotime values, where convergence is defined as correlation of cell indices along the trajectory.

I abandoned the project after observing that more iterations to not change the convergence towards the actual solution. While pseudotime values close to 0 recapitulated the actual solution quite well, pseudotime values towards 100 failed to correlate at all with the actual solution. There was no convergence with more iterations.

Maybe this will be helpful to somebody, which is why I'm posting.

# Parallel Iterative Sampling with Slingshot
# pissshot is a course-grained estimator of the actual slingshot solution, and iteratively converges towards the actual solution
# pissshot randomly partitions cells in a large database into bins, bin size is determined by finding the computationally optimal number of cells
# pissshot is run independently on each bin, pseudotime vectors are aligned across all bins, and then all bins are merged
# pissshot repeats this process and the average of pseudotime assignments for each cell is taken at the end of each iteration.
# The convergence of the blingshot model at the end of each iteration is calculated by running slingshot on a small sample of adjacent cells at a random pseudotime interval and measuring correlation with the predicted model
# When satisfactory convergence is reached, blingshot returns a singlecellobject, just as slingshot would

pissshot <- function(sce, num_iterations = 20){
    cat("Running blingshot on", dim(sce)[2],"cells and",dim(sce)[1],"features\n")
    cat("  Step 1/5: Finding computationally optimal increment sizes\n   ")
#    increment_size <- FindOptimalIncrementSize(sce)
    increment_size <- 688
    num_increments <- ceiling(dim(sce)[2]/increment_size)-1
    cat("\n  ...Increment size of",increment_size,"cells is computationally optimal\n")
    cat("  Step 2/5: Course-grained pseudotime assignment across",num_increments,"increments")
    ps <- list()
    convergenceArr <- c()
    for(iteration in seq(from=1,to=num_iterations,by=1)){
        cat("\n   ... iteration",iteration,"...\n")
        cat("         randomly assigning cells to",num_increments,"bins\n         0%.")
        bins <- PartitionBins(sce, increment_size = increment_size, seed = 7*iteration)
        cat("100%\n         running slingshot on each bin\n         0%")
        cell_pseudotimes <- c()
        for(b in seq(from=1,to=length(bins),by=1)){
            bin_with_pseudotime <- suppressMessages(slingshot(bins[[b]]))
            # Find right cell and left cell in pseudotime array (cell with highest/lowest pseudotime values)
            pcells <- as.data.frame(colData(bin_with_pseudotime)$slingPseudotime_1)
            rownames(pcells) <- colnames(bin_with_pseudotime)
            pcells <- data.frame(lapply(pcells, function(x) as.numeric(as.character(x))), check.names = F, row.names = rownames(pcells))
            pcells.sorted <- pcells[order(-pcells[,1]), , drop = FALSE]
            right_cell_ID <- rownames(pcells.sorted)[dim(pcells.sorted)[1]]
            left_cell_ID <- rownames(pcells.sorted)[1]
            right_cell_pos <- match(right_cell_ID,colnames(bins[[b]]))
            left_cell_pos <- match(left_cell_ID,colnames(bins[[b]]))
            right_cell_logcounts <- assays(bins[[b]])$logcounts[,right_cell_pos]
            left_cell_logcounts <- assays(bins[[b]])$logcounts[,left_cell_pos]
            if(b==1&&iteration==1){
                ref_right_cell_logcounts <- right_cell_logcounts
                ref_left_cell_logcounts <- left_cell_logcounts
            }
            # Figure out whether this pseudotime in this bin has flipped relative to reference (first bin), and if so, flip pseudotime values
            cor_opposite_ends <- mean(cor(right_cell_logcounts,ref_left_cell_logcounts),cor(left_cell_logcounts,ref_right_cell_logcounts))
            cor_same_ends <- mean(cor(right_cell_logcounts,ref_right_cell_logcounts),cor(left_cell_logcounts,ref_left_cell_logcounts))
            flipped <- FALSE
            if(cor_opposite_ends > cor_same_ends){
                flipped <- TRUE
            }
            if(flipped==TRUE){
                # find max pseudotime in colData(bin_with_pseudotime)$slingPseudotime
                # Recalculate pseudotime in bin_with_pseudotime as max-value
                pseudotimes <- colData(bin_with_pseudotime)$slingPseudotime_1
                colData(bin_with_pseudotime)$slingPseudotime_1 <- max(pseudotimes)-pseudotimes
            }
            cells_with_pseudotime <- as.data.frame(colData(bin_with_pseudotime)$slingPseudotime_1,colnames(bin_with_pseudotime))
            cell_pseudotimes <- rbind(cell_pseudotimes, cells_with_pseudotime)
            cat(".")
        }
        cat("100%\n")

        ps[[iteration]] <- cell_pseudotimes[order(row.names(cell_pseudotimes)), , drop = FALSE]

        pss <- data.frame(ps)
        if(iteration>1){
            convergence10 <- measureConvergence(pss, sce = sce, ptime = 10, increment_size = increment_size)
            convergence50 <- measureConvergence(pss, sce = sce, ptime = 50, increment_size = increment_size)
            convergence99 <- measureConvergence(pss, sce = sce, ptime = 99, increment_size = increment_size)
            cat("   convergence at 10:",convergence10,"... at 50:",convergence50,"...at 100:",convergence99,"\n")
            convergenceArr <- rbind(convergenceArr,c(iteration,convergence10,convergence50,convergence99))
        }
        # Iterate until the model converges towards the actual slingshot solution, within an indicated fraction (i.e. convergence = 0.02)
        # The convergence measure compares how similar the model is to the actual slingshot solution, i.e. a convergence of 0.02 means that cell pseudotime values are within 2% of the actual solution
        # Assess this by selecting a pseudotime at random, analyzing the following increment_size cells, and comparing accuracy of pseudotime measurement across that trajectory with the composite of the course-grained analysis

        # need to study how convergence changes over pseudotime
    }

    # remove cells with really high standard deviations

    return(pss)
#    return(convergenceArr)
}

measureConvergence <- function(mat, sce = sce, ptime = 45, increment_size = 500){
    # this line is throwing an error
    mat[,"avg"] <- apply(mat[,1:dim(mat)[2]],1,mean)
    mat <- mat[order(mat[,"avg"]), , drop = FALSE]
    # select a random pseudotime increment
    mat <- mat[mat[, "avg"] > ptime,]
    mat <- mat[1:increment_size,]
    cell_IDs <- rownames(mat)
    # pull out this list of cell_IDs from the sce object and run slingshot
    sce2 <- suppressMessages(slingshot(sce[,cell_IDs]))
    cpcells <- as.data.frame(colData(sce2)$slingPseudotime_1)
    rownames(cpcells) <- colnames(sce2)
    cpcells <- data.frame(lapply(cpcells, function(x) as.numeric(as.character(x))), check.names = F, row.names = rownames(cpcells))
    # sort cpcells by increasing pseudotime
    cpcells <- cpcells[order(cpcells[,1]), , drop = FALSE]
    ordered_cells <- cbind(rownames(mat),rownames(cpcells))
    # get the index of each cell in the dataframe, run a correlation on how well the indices line up
    cell_indices <- c()
    for(cell in cell_IDs){
        new_indices <- c(which(ordered_cells[,1] == cell), which(ordered_cells[,2] == cell))
        cell_indices <- rbind(cell_indices,new_indices)
    }
    convergence_measure <- cor(cell_indices[,1],cell_indices[,2])
    plot(cell_indices[,1],cell_indices[,2])
    return(convergence_measure)
}

PartitionBins <- function(sce, increment_size = 500, seed = 123){
    bins <- c()
    while(dim(sce)[2]>=increment_size){
        if(dim(sce)[2]<increment_size*2){
            bins <- c(bins, sce)
            sce <- sce[,1]
        }
        else{
            set.seed(seed)
            train_ind <- sample(seq_len(ncol(sce)), size = increment_size)
            bins <- c(bins, sce[,train_ind])
            sce <- sce[,-train_ind]
            cat(".")
        }
    }
    return(bins)
}

FindOptimalIncrementSize <- function(sce){
    # Find the computationally optimal increment_size for grainy alignments using 5 increment points (100, 500, 1000, 2500)
    res <- c(100,250,500,750,1000,1500,2500,5000)
    runtimes <- c()
    for(i in res){
        runtimes <- c(runtimes, suppressMessages(system.time(slingshot(sce[,0:i])))["elapsed"])
    }
    mat <- as.data.frame(t(rbind(res,runtimes)))
    c <- summary(lm(mat$runtimes ~ poly(mat$res, 2, raw=TRUE)))$coefficients[,"Estimate"]
    predicted_runtimes <- c()
    for(i in seq(from=100, to=5000, by=1)){
        predicted_runtime <- (c[1] + c[2]*i + c[3]*i^2)*(dim(sce)[2]/i)
        predicted_runtimes <- rbind(predicted_runtimes,c(i,predicted_runtime))
        cat(".")
    }
    increment_size <- predicted_runtimes[which.min(predicted_runtimes[,2]),1]
    return(increment_size)
}


After completing the evaluation, how to predict PTM lysine or any other modification with unlabeled data for making a web server? I have created a model with SVM classifier and evaluated it with 5-fold cross-validation. The dataset was sequence data containing lysine modified sites. Now I want to build a web server where anyone can input any sequence in FASTA format. Apparently, I have no label information at this point. So how can I will be able to extract features or use any feature selection method which requires label information? How can I extract or create the test set to predict? Suppose I have used f-score feature selection method in my model and selected 50 features from the given dataset. But for web server, having no label information, how can I select 50 features from the sequence that are provided by an user and then test it with my previously trained model?

As being a novice in the world of bioinformatics, I am seeking the answer which might be improper for others. But it will be very much helpful for me.
I am giving some link of some papers and its corresponding web servers.

iPTM-mLys: identifying multiple lysine PTM sites and their different type (https://academic.oup.com/bioinformatics/article/32/20/3116/2196525)

Web server link : http://www.jci-bioinfo.cn/iPTM-mLys
iHyd-PseCp: Identify hydroxyproline and hydroxylysine in proteins by incorporating sequence-coupled effects into general PseAAC
(https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5190098/)

Web server link : http://www.jci-bioinfo.cn/iHyd-PseCp

Interesting application. I would suggest using Python's Flask system, which is described in a few pages in O'Reilly's Programming PyTorch, chapter 8. PyTorch in Production. Django is not really recommended (don't say that on Stackoverflow however). O'Reilly do an entire book on Flask, but thats to make the site pretty.

Obviously, the explanation here is Deep Learning rather than ML but the basic gist will hold. You should also think about building a docker for distribution which is easy (same chapter). We're in lockdown so I don't have this book with me.

The difficult part will be web architecture because you shouldn't have the CPU/GPU on the same server as the quiry interface, otherwise you are prone to a denial of service attack or just sytem overload. To do this you need to consider  Representational state transfer (RESTful) and that is not easy at all (at least for me). Flask facilitates RESTful, but thats like systems admin stuff. I guess before getting into RESTful I'd see how popular your website is, but if its the same server as your institute website you'd need to be really careful because an external user could bring your institutional website down (you wouldn't be popular).


Speeding up network creation in WGCNA I am following the tutorials on the WGCNA website for creating a coexpression network. It is taking way too long. Is there a way to speed up execution? All I need is to create the co-expression graph. I will then export to Python and do my analysis there.
In particular, the step "calculating module eigengenes" is taking a lot of time. I do not really care about the modules. All I want is the graph with the genes and the weighted edges.

You can use the "manual network construction" approach. Calculate TOM, then use exportNetworkToCytoscape or exportNetworkToVisANT to export the network (or simply save the TOM matrix as a table and load it in Python).


Hierarchical models with limma? I have a dataset with (microarray) gene expression data that was sampled from the same individuals at multiple timepoints. Our exposure is a continuous variable, and because this was an observational study there is no consistent pattern of change in exposure over time (eg. one person's exposure may continuously increase, another's may go up and down, another may have no changes, etc.). The exposure is not expected to have any cumulative effect, so different timepoints within a person do not influence each-other. Our goal is to find genes where expression for the gene is associated with the level of exposure.

We are thinking of analyzing this as a hierarchical model, with timepoints nested within person (so a mixed model with a random effect for person). We will also need to control for a few covariates, some of which will be the same at all timepoints (race, sex) and some of which may change between timepoints (use of particular medications).

Is this something that it is possible to use limma for? I found that limma has a function called duplicateCorrelation() intended for block designs that some people have tried to use to fit mixed models, but I am not sure whether it can be used for this. Or if not limma, then what would you suggest?

Yes, you can use limma for this mixed model approach. Like you suggest, the random effect (persons) can be put in duplicateCorrelation().

Here is a similar example with RNAseq data, on bioconductor support site.


Why aren't clustal alignments stable under deletion of some of the sequences? I'm new to pairwise and multiple sequence alignment in general, but I thought I understood how clustal works -- k-tuple distance is a cheap metric that's 'probably approximately' monotonic in the real global pairwise-alignment score, so we can just calculate it between all C(n, 2) sequences and do a simple UPGMA. Clustalw algorithm description here

But I'm observing odd results. Suppose I have a file with 60 sequences and run clustal on it, and observe that seqA and seqB are paired terminal child nodes. If I then delete a handful (say 10) of the sequences, but keep seqA and seqB, then in the new clustal alignment, seqA and seqB may not remain paired.

This is strange to me, because if Clustal was actually doing the k-tuple distance + UPGMA thing I thought it was doing, then all terminal pairings should always be stable under deletions of other sequences.

So... what am I missing?

As the graph of sequences depends on the sequences in the input, if you change the input sequences (by deleting some sequences) you change the graph.
The alignment between the sequences changes, and thus the relative distance between them changes too.

The deletion of some sequences won't cause any change on other sequences only if the full space of the C(n 2) pairwise comparisons would be done. However for a large number (and long) of sequences, this would require too much computational effort and time, making such methods as the k-tuple distance competitive.


A new paper suggests the Corona Virus has "Uncanny similarity of unique inserts in the 2019-nCoV spike protein to HIV-1" - What does this mean? Quote:

  We found 4 insertions in the spike glycoprotein (S) which are unique to the 2019-nCoV and are not present in other coronaviruses. Importantly, amino acid residues in all the 4 inserts have identity or similarity to those in the HIV-1 gp120 or HIV-1 Gag. Interestingly, despite the inserts being discontinuous on the primary amino acid sequence, 3D-modelling of the 2019-nCoV suggests that they converge to constitute the receptor binding site. The finding of 4 unique inserts in the 2019-nCoV, all of which have identity/similarity to amino acid residues in key structural proteins of HIV-1 is unlikely to be fortuitous in nature.

Source: https://www.biorxiv.org/content/10.1101/2020.01.30.927871v1

I know this is a preprint paper and it's not yet peer reviewed, but can someone tell me what the implications are if this is true? Does this mean that the virus is artificially constructed?

UPDATE: The article has now been withdrawn with the following note:

  This paper has been withdrawn by its authors. They intend to revise it
  in response to comments received from the research community on their
  technical approach and their interpretation of the results. If you
  have any questions, please contact the corresponding author.

This is very odd, and will require a rigorous investigation, but my initial reaction is one of scepticism. Considering just the 1st insert, the insert sequence is GTNGTKR, short at just 7 amino acids. A simple BLASTP vs NR did not find perfect matches to HIV sequences, but did reveal 100% identity across the full sequence length to >50 other short protein sequences, which could be spurious 'chance' hits of course. Many of the top 100 hits are against eukaryotic sequences, for example one is against Pristionchus pacificus, a type of nematode worm (see alignment below). Someone needs to do a proper peer review of this preprint before any conclusions are drawn from it.

>tank-1 [Pristionchus pacificus]
Sequence ID: PDM74036.1 Length: 2481
Range 1: 1474 to 1480

Score:24.0 bits(49), Expect:2477,
Method:,
Identities:7/7(100%), Positives:7/7(100%), Gaps:0/7(0%)

Query  1     GTNGTKR  7
             GTNGTKR
Sbjct  1474  GTNGTKR  1480

UPDATE 1: Things are moving fast. There are now 10 comments under the preprint, all of them critical of the idea that these inserts are meaningfully similar to HIV-1 sequences for the reason I outlined above (i.e. short sequences with many hits against organisms from across the tree of life) and additional points, such as small insertions being quite normal evolution for RNA viruses. There are also similar critiques on twitter and biorxiv has added the following disclaimer header to their website:

  bioRxiv is receiving many new papers on coronavirus 2019-nCoV.   A
  reminder: these are preliminary reports that have not been
  peer-reviewed. They should not be regarded as conclusive, guide
  clinical practice/health-related behavior, or be reported in news
  media as established information.

I have never seen such rapid (and unanimous) post-publication peer review before. It is therefore already (just 1 day after publication) fairly clear that at least some claims in this study are completely false and there is no good evidence to support the claim that the 2019-nCoV have acquired sequences from HIV or evidence to suggest the virus is engineered.


How to obtain NCBI Nucleotide (Entrez) IDs for all E coli sequences? I would like to gather a list of ALL E coli NCBI (Entrez/Nucleotide) IDs so that I can download fna and faa files, through a function like this:
# python code #
from Bio import Entrez
from glob import glob

def dl_genome(id, folder='genomes'):
    files=glob('%s/*.fna'%folder)
    out_file = '%s/%s.fna'%(folder, id)

    if out_file in files:
        print (out_file, 'already downloaded')
        return
    else:
        print ('downloading %s from NCBI'%id)

    Entrez.email = ""
    handle = Entrez.efetch(db="nucleotide", id=id, rettype="fna", retmode="text")
    with open(out_file, 'w') as fout:
        fout.write(handle.read())

for strain in targetStrainIDs:
    dl_genome(strain, folder='genomes')

My problem is that I don't have a good way to get the E coli strain names and their corresponding NCBI IDs for Entrez. Is there a way to get these NCBI IDs without manually searching for them?
Context
I have been trying to develop an E coli pan-genome for one of my projects. I initially just downloaded the relevant genomes (fna, faa, & features.tab files) from the PATRIC database. However, I've come to realize that the reference sequence on PATRIC has issues. Specifically, the reference sequence on PATRIC's fna/faa files for E coli K12 MG1655 notably do NOT list genes like b2031 or b2033, despite both being present in this strain.
The NCBI reference genome DOES include these genes (with proper annotations), so I'd like to just download what I need from NCBI instead. However, I'd like a systematic way to download all the fna/faa files since my pan-genome is eventually going to be built from hundreds of of these fna/faa files.
Any help would be appreciated.

Do you mean the NCBI taxonomy IDs? If so, I've personally always preferred the way UniProt displays NCBI taxonomy, and it's easier to download IDs in mass.
For example, if you find E.coli's entry: https://www.uniprot.org/taxonomy/562
You can then click on "all lower taxonomy nodes", which will show all E.coli strains. If you click 'Preview first 10' before downloading, you'll see that the first column is the NCBI taxonomy ID of each strain of E.coli.
For example, E.coli strain 0.1197 is reported to have the ID 1005443 from this preview page. Indeed, if we search this ID on NCBI's taxonomy page, we do get the correct entry: https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=1005443
If this is what you want, click 'download all' to get a nice tab-separated file of all the strains of E.coli and their NCBI taxonomy IDs. This file could then be read and processed before you send API requests to Entrez.
Hope this helps!


How can I interpret gene expression data from Bioconductor packages? I am currently looking at microarray data from a bioconductor microarray dataset. Specifically, I have data (a snippet) which looks like the following:

   GSM794885   GSM794886   GSM794887   GSM794888   GSM794889
   -1.67170290 -0.07803822  0.99990560 -1.94692850 -1.78671360

I am wondering how to interpret the values in this microarray dataset. Is it generally the case that larger values mean a more significant effect?

Processed microarray data is commonly represented as log_2(expression). This data transformation is used because the data more closely fit a normal distribution in log space. With such a transformation, negative values represent very small amounts, rather than the absence of something.

Larger values mean that more stuff is being produced (presumably transcripts), but its significance can only be assessed when combined with other data (e.g. through differential expression analysis).


Calculating the charge of a peptide computationally I was wondering how I can calculate the charge of a protein peptide (e.g. "RKTTLVPNTQTASPR") computationally in R or another tool.

A quick google search turns up protcalc, which is able to give a nice pH-dependent table of peptide charges (yours ranges from 3.1 at pH 4 to 1.5 at pH 10). It's on sourceforge, so hopefully the source code (and maybe a publication) is somewhere in there. Granted, this isn't R, but it's not clear from your post how necessary that really is.

BTW, in general "best and most accurate" turns out to be highly subjective (at least the "best" part of that).


How can I calculate the silent mutation for each position in codon? I am trying to find the number of silent mutations in for each position in codon.

In the homework, I am asked:

  For the codons and for each position, count the number of silent mutation. There are 64 codons, for each position there are 3 mutations leading to a change at the DNA level. There are thus 192 possible mutations for each of the three positions of the codons. For each position, print the number of silent mutations.

I understand stand silent mutation is the mutation that doesn't lead to amino acid change. For example, GCC can become GCT but it is still Thr

I am not sure what does "calculating silent mutation for each position" mean? I mean it is easy to find silent mutation for each codon like the example I gave above, but couldn't really understand how to calculate the silent mutation for each position.

Should I change letters (nucleotides) one by one and see if it causing a silent mutation?

Calculate in this context would just mean determine. You don't actually have to change every base to determine whether a change is silent. Instead, try building a hash table of amino acids to codons, which you can then quickly iterate over to find all silent changes.

BTW, try to avoid loaded terms like "mutation", even if your professor incorrectly used it. Terming things like this as "silent polymorphism" or "silent change" or something along those lines is preferable.


Percentage distribution of cells in all clusters based on their treatment condition? I have 2151 cells, I clustered them by Seurat to 5 clusters. With the code below, I am able to have the number of cells per cluster and per condition:

number_perCluster<- table(object@meta.data$conditions,
                          object@meta.data$clusterID)

Could I please get a hint on how I could proceed to produce a figure where I say: X% of cells treated with condition "Y" are located in cluster "C" and so forth?

Thank you in advance.

Here is a solution using dplyr and ggplot2:
library(Seurat)
library(dplyr)
library(ggplot2)

meta.data <- pbmc_small[[]]

# create random classifications for the sake of this example
meta.data$condition <- sample(c('A', 'B', 'C'), nrow(meta.data), replace = TRUE)

counts <- group_by(meta.data, condition, res.1) %>% summarise(count = n())

ggplot(counts, aes(res.1, count, fill = condition)) +
  geom_bar(stat = 'identity')


Visualising gene expression across cell type and conditions in one plot, in Single Cell Sequencing data I am new to single cell sequencing data analysis but I have basic programming background in R and python. I want to be able to plot differential expression of two genes, Gene1 and Gene2, across three cell types and across three conditions on one plot. I have followed this tutorial to perform QC and various differential analyses using Seurat package on my data up till now.

My data set is one object, containing information about the three cell types under three various conditions. The metadata slot of my data set contains information about my cell types as well as the conditions under which they are tested. Using the following DotPlot commands I am able to generate separate plots of gene expression with respect to cell type and with respect to condition:

Seurat::DotPlot(object = object, genes.plot = c("Gene1", "Gene2"), plot.legend = TRUE, group.by = "cell_types")

Seurat::DotPlot(object = object,
        genes.plot = c("Gene1", "Gene2"), plot.legend = TRUE, group.by = "tissue_types" )

For elegant representation of my data, I want to be able to have the information of the following plots into one plot.

I will be very grateful for any hint on how I can proceed to achieve an elegant presentation as well as information about other possible ways to graphically represent gene expression data for publication or presentation purposes.

Thank you in advance for your help.

Try using SplitDotPlotGG, this will allow you to color the dots by some additional variable (eg condition).

In your case you would want to do something like this:

library(Seurat)

SplitDotPlotGG(object = object,
               genes.plot = c("Gene1", "Gene2"),
               plot.legend = TRUE,
               group.by = "cell_types",
               grouping.var = "tissue_types")

This would give a result similar to below (an example from https://satijalab.org/seurat/immune_alignment.html).


RNA-Seq data transformation prior to sample correlation analysis If I’m starting with Deseq2 normalized counts what are some preprocessing steps that I should apply to these data before estimating sample correlation using the cor function in R? For example, would it make sense to quantile normalize the data?

You don't absolutely need to perform any further steps before calculating the correlation. However you would likely benefit from the following sort of workflow:

Transform your data with vst() or robust log (see the DESeq2 vignette) so genes are on a more or less similar scale.
Select the top 300-500 most variable genes
Use this subset with to compute sample-wise correlation.

This also what follows what you would find in the DESeq2 user guide, which I strongly encourage you to read.


BED file from .bam alignment structure I am using a pipeline from nanopore to detect structural variants (SVs) in a human sample with long-reads sequencing. The first steps of the pipeline are:

index the reference genome with minimap2
map the long reads to a reference genome using same tool as above
Get BED coordinates from the .bam file.

While the first two steps uses minimap, the third one uses a python script written by the developers of the pipeline, which can be found here.

After the pipeline calls bamref2bed i.e. the python script above, my target.bed looks like the following:

chr1    0   249250621
chr2    0   243199373
chr3    0   198022430
chr4    0   191154276
chr5    0   180915260
chr6    0   171115067
chr7    0   159138663
chr8    0   146364022
chr9    0   141213431
chr10   0   135534747
chr11   0   135006516
chr12   0   133851895
chr13   0   115169878
chr14   0   107349540
chr15   0   102531392
chr16   0   90354753
chr17   0   81195210
chr18   0   78077248
chr19   0   59128983
chr20   0   63025520
chr21   0   48129895
chr22   0   51304566
chrX    0   155270560
chrY    0   59373566

As far as I am concerned the first three columns of a bed file refer to (I) chr name, (II) start position (III) end position. The main thing I am wondering is how come all chromosomes have a start position of 0. Also, reading about converting bam to bed files I guess the most important feature would be the coverage across the whole chromosomes, which in this file is missing.

But my main question would be: is it normal to have a start position of 0 for all chromosomes or is there something wrong here?

I cannot access the pipeline from nanopore, but looking at the python script you provided:

with pysam.AlignmentFile(bam_file, "r") as bam:
    for tid in range(0, bam.nreferences):
        ref_name = bam.get_reference_name(tid)
        if len(list(filter(lambda x: x in ref_name, filter_names))) > 0:
            continue

        print(ref_name,
              "0",
              str(bam.get_reference_length(ref_name)),
              sep='\t')

So the script gets the length of the chromosome and by default, outputs the start as 0.Zero based coordinate system is sometimes used in UCSC systems, so I am only guessing the authors of the pipeline use this convention somewhere further down.

As for the coverage you mentioned, that is usually output as a score in the 5th column (bed format). In those situations, usually the coordinates of the alignments in the bam file are converted into bed format. In this pipeline, i think the authors are simply extracting the chromosome information out and putting it into a tabular form.


Is removing samples based on clustering for downstream analysis a right choice? I'm using TCGA Lung cancer data. I'm interested in doing differential analysis between Lung vs Normal. Before DEA, to check the distance between each pair of samples I plotted an MDS plot:

In this, I see some Tumor samples are clustered with Normals. As this is a two dimensional plot, it is a little hard to say.

So, I made a clustering heatmap with the top 10% genes:

 (Where brown color is Normal and black color is Tumor). I also see some Normals clustering with Tumor here, and some Tumor are with the Normals.

My question is: Should I exclude those tumor samples which are clustering with the normals for differential analysis? Is it a good idea to remove samples based on clustering?

Note: I didn't see any technical bias of those samples

I would be very hesitant to blindly exclude those samples based on the clustering. Check to see if the clusters actually denote some sort of batch effect, since it's not like all of the TCGA datasets were processed at the same place or time. Check for clinical covariates too. Also, do a PCA and check the genes with high loading on the relevant PC. If it appears that the genes are indicators for cancerous state, then it's likely the samples are just mislabeled and you can reasonably exclude them. Otherwise keep them in so your results better represent the underlying biological variability.

BTW, your heatmap and MDS plots are what I'd expect from a cancer dataset. You occasionally see a tumour with the normal samples and vice versa, since the problematic genes aren't always in the top 10% and there is likely some heterogeneity in the cancer samples.


Regressing out unwanted sources of variation in single cell RNA-seq data I have a dataset of single cell count data and I want to regress out the variation caused by the number of UMI's and the percentage of mitochondrial genes.

I know that count data is discrete data, and commonly follows a negative binomial distribution. To regress out these two confounders, should I use a linear model or GLM (poisson/negative binomial)? And how can I determine the optimal choice?

Cleaning data before doing analysis can be more important than the method of data analysis. For UMIs, there's an obvious cleaning step that can be done: it would be better to filter out duplicates (prior to generating counts) than to try to incorporate that information in the count analysis.

I'm not familiar with single cell analysis, but I am aware that DESeq2 can be adapted to work with single cells, and incorporating additional variation is just a matter of adding an additional covariate into the DESeq2 formula model.


How can I create my own GO association file (gaf)? This question is based on a question on BioStars  posted >2 years ago by user jack.

It describes a very frequent problem of generating GO annotations for non-model organisms. While it is based on some specific format and single application (Ontologizer), it would be useful to have a general description of the pathway to getting to a GAF file.

Note, that the input format is lacking a bit of essential information, like how it was obtained. Therefore, it is har to assign evidence code. Therefore, lets assume that the assignments of GO terms were done automagically.

  I want to do the Gene enrichment using Ontologizer without a
  predefined association file(it's not model organism).

  I have parsed a file with two columns for that organism like this :

geneA  GO:0006950,GO:0005737
geneB  GO:0016020,GO:0005524,GO:0006468,GO:0005737,GO:0004674,GO:0006914,GO:0016021,GO:0015031
geneC  GO:0003779,GO:0006941,GO:0005524,GO:0003774,GO:0005516,GO:0005737,GO:0005863
geneD  GO:0005634,GO:0003677,GO:0030154,GO:0006350,GO:0006355,GO:0007275,GO:0030528

  I have downloaded the .ob file from Gene ontology file which contain
  this information (from here) :

!
! GO IDs (primary only) and name text strings
! GO:0000000 [tab] text string [tab] F|P|C
! where F = molecular function, P = biological process, C = cellular component
!
GO:0000001  mitochondrion inheritance   P
GO:0000002  mitochondrial genome maintenance    P
GO:0000003  reproduction    P
GO:0000005  ribosomal chaperone activity    F
GO:0000006  high affinity zinc uptake transmembrane transporter activity    F
GO:0000007  low-affinity zinc ion transmembrane transporter activity    F
GO:0000008  thioredoxin F
GO:0000009  alpha-1,6-mannosyltransferase activity  F
GO:0000010  trans-hexaprenyltranstransferase activity   F
GO:0000011  vacuole inheritance P

  What I need as output is .gaf file in the following format (in the
  format of the files here):

!gaf-version: 2.0

!Project_name: Leishmania major GeneDB

!URL: http://www.genedb.org/leish

!Contact Email: mb4@sanger.ac.uk

 GeneDB_Lmajor    LmjF.36.4770    LmjF.36.4770        GO:0003723    PMID:22396527    ISO    GeneDB:Tb927.10.10130    F    mitochondrial RNA binding complex 1 subunit, putative    LmjF36.4770    gene    taxon:347515    20120910    GeneDB_Lmajor
 GeneDB_Lmajor    LmjF.36.4770    LmjF.36.4770        GO:0044429    PMID:20660476    ISS        C    mitochondrial RNA binding complex 1 subunit, putative    LmjF36.4770    gene    taxon:347515    20100803 GeneDB_Lmajor             GeneDB_Lmajor    LmjF.36.4770    LmjF.36.4770        GO:0016554    PMID:22396527    ISO    GeneDB:Tb927.10.10130    P    mitochondrial RNA binding complex 1 subunit, putative    LmjF36.4770    gene   taxon:347515    20120910    GeneDB_Lmajor
 GeneDB_Lmajor    LmjF.36.4770    LmjF.36.4770        GO:0048255    PMID:22396527    ISO    GeneDB:Tb927.10.10130    P    mitochondrial RNA binding complex 1 subunit, putative    LmjF36.4770    gene    taxon:347515    20120910    GeneDB_Lmajor

  How to create your own GO association file (gaf)?

Here's a Perl script that can do this:

#!/usr/bin/env perl
use strict;
use warnings;

## Change this to whatever taxon you are working with
my $taxon = 'taxon:1000';
chomp(my $date = `date +%Y%M%d`);

my (%aspect, %gos);
## Read the GO.terms_and_ids file to get the aspect (sub ontology)
## of each GO term.
open(my $fh, $ARGV[0]) or die "Need a GO.terms_and_ids file as 1st arg: $!\n";
while (<$fh>) {
    next if /^!/;
    chomp;
    my @fields = split(/\t/);
    ## $aspect{GO:0000001} = 'P'
    $aspect{$fields[0]} = $fields[2];
}
close($fh);

## Read the list of gene annotations
open($fh, $ARGV[1]) or die "Need a list of gene annotattions as 2nd arg: $!\n";
while (<$fh>) {
    chomp;
    my ($gene, @terms) = split(/[\s,]+/);
    ## $gos{geneA} = (go1, go2 ... goN)
    $gos{$gene} = [ @terms ];
}
close($fh);

foreach my $gene (keys(%gos)) {
    foreach my $term (@{$gos{$gene}}) {
        ## Warn and skip if there is no aspect for this term
        if (!$aspect{$term}) {
            print STDERR "Unknown GO term ($term) for gene $gene\n";
            next;
        }
        ## Build a pseudo GAF line
        my @out = ('DB', $gene, $gene, ' ', $term, 'PMID:foo', 'TAS', ' ', $aspect{$term},
                             $gene, ' ', 'protein', $taxon, $date, 'DB', ' ', ' ');
        print join("\t", @out). "\n";
    }
}

Make it executable and run it with the GO.terms_and_ids file as the 1st argument and the list of gene annotations as the second. Using the current GO.terms_and_ids and the example annotations in the question, I get:

$ foo.pl GO.terms_and_ids file.gos
DB  geneD   geneD       GO:0005634  PMID:foo    TAS     C   geneD       protein taxon:1000  20170308    DB
DB  geneD   geneD       GO:0003677  PMID:foo    TAS     F   geneD       protein taxon:1000  20170308    DB
DB  geneD   geneD       GO:0030154  PMID:foo    TAS     P   geneD       protein taxon:1000  20170308    DB
Unknown GO term (GO:0006350) for gene geneD
DB  geneD   geneD       GO:0006355  PMID:foo    TAS     P   geneD       protein taxon:1000  20170308    DB
DB  geneD   geneD       GO:0007275  PMID:foo    TAS     P   geneD       protein taxon:1000  20170308    DB
DB  geneD   geneD       GO:0030528  PMID:foo    TAS     F   geneD       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0016020  PMID:foo    TAS     C   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0005524  PMID:foo    TAS     F   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0006468  PMID:foo    TAS     P   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0005737  PMID:foo    TAS     C   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0004674  PMID:foo    TAS     F   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0006914  PMID:foo    TAS     P   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0016021  PMID:foo    TAS     C   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0015031  PMID:foo    TAS     P   geneB       protein taxon:1000  20170308    DB
DB  geneA   geneA       GO:0006950  PMID:foo    TAS     P   geneA       protein taxon:1000  20170308    DB
DB  geneA   geneA       GO:0005737  PMID:foo    TAS     C   geneA       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0003779  PMID:foo    TAS     F   geneC       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0006941  PMID:foo    TAS     P   geneC       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0005524  PMID:foo    TAS     F   geneC       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0003774  PMID:foo    TAS     F   geneC       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0005516  PMID:foo    TAS     F   geneC       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0005737  PMID:foo    TAS     C   geneC       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0005863  PMID:foo    TAS     C   geneC       protein taxon:1000  20170308    DB

Note that this is very much a pseudo-GAF file since most of the fields apart from the gene name, GO term and sub-ontology are fake. It should still work for what you need, however.


GFF file too big to load into R I want to load a GFF (annotation) file into an R data frame in order to extract some information from it, e.g. the location of transposons in the genome.

I used the following code:

library(rtracklayer)
gff <- rtracklayer::import('c_elegans.PRJNA13758.WS263.annotations.gff3.gz')
Error: cannot allocate vector of size 210.2 Mb

It seems that the GFF3 file is too big for loading into R. Do you know if this is limited by R or by my computer's memory?

I ran this on a GNU/Linux computer.

It depends of course on your available RAM on your computer.

Under Windows you can increase the available RAM in R with memory limit:

> memory.limit()
[1] 8070
> memory.limit(size=10000)
[1] 10000
> memory.limit()
[1] 10000

But that is not necessary on linux/unix.


EGA remote download via snakemake fails due to EGA_CLIENT_ID specification The example code provided by the Snakemake documentation for remote downloading from the EGA archive fails for me. I have defined all required environment variables, but an error comes up stating

{"error":"invalid_client","error_description":"Bad client credentials"}

I guess I am not defining the EGA_CLIENT_ID variable correctly? I am using Snakemake 5.16 and Python 3.7.2. Thanks a lot in advance for your advice.

Code:

import snakemake.remote.EGA as EGA

ega = EGA.RemoteProvider()
os.environ["EGA_USERNAME"] = "XXX"
os.environ["EGA_PASSWORD"] = "XXX"
os.environ["EGA_CLIENT_SECRET"] = "XXX"
os.environ["EGA_CLIENT_ID"] = "pyega3"
rule a:
  input:
      ega.remote("ega/EGAD00001002142/COLO_829_EPleasance_TGENPipe.bam.bai")
  output:
      "data/COLO_829BL_BCGSC_IlluminaPipe.bam.bai"
  shell:
      "cp {input} {output}"

Okay, I finally solved it. The reason was the EGA_CLIENT_ID variable. In the pyega3.py script, there is a definition for the client_id for pyega3 v3.0.44 ("f20cd2d3-682a-4568-a53e-4262ef54c8f4"). After setting the EGA_CLIENT_ID to this value login was successful.

os.environ["EGA_CLIENT_ID"] = "f20cd2d3-682a-4568-a53e-4262ef54c8f4"

I think this is a bug and I filed it on Snakemake's git repository - the official documentation just requires definition of the EGA_USERNAME and EGA_PASSWORD variables, and while definition of EGA_CLIENT_SECRET will probably need to be done by the user, the client ID should be pulled automatically from the pyega3 version installed locally.

Of note, after successful login, the example download from the documentation failed with the following:

Missing input files for rule a:
ega/EGAD00001002142/COLO_829_EPleasance_TGENPipe.bam.bai

However, downloading files that I have access to worked, indicating that the paths to the files from the documentation are wrong or are access-protected.


Identify differentially covered genes only between two samples I have a question about finding differentially covered regions (coverage represents methylation level which goes from 0 to several thousands). I'm using enrichment based method which can be summarized with coverage per gene:

data <- matrix(sample(80), 20)
# Genes in rows
rownames(data) <- letters[1:20]
colnames(data) <- c("group_A_tr1", "group_A_tr2", "group_B_tr1", "group_B_tr2")

In data matrix each row represents a gene and each column represents a sample. There are two sample groups (A and B) with two technical replicates per each group. Problem is that we do not have any biological replicates.

My goal is to identify genes that are differentially methylated between two groups. I know that limma, edgeR, DESeq2 can be used in analysis like this, however I don't have enough samples. Basically I'll need to compare only two columns (after averaging technical replicates).

What method would be appropriate to work with data like this? Is it possible to treat technical replicates as biological ones?

Not really an answer but an extended comment... and most likely something you don't like to hear

I guess by technical replicates, it means taking the same biological sample and making 2 methylation libraries. If this is used as replicates in deseq2 or edger, the variance you are estimating is the technical variation that comes with preparing the library, bisulphite conversion etc..

Often the biological variance will be much larger. That is if you take 2 independent samples of group A, their variation in methylation is much larger than the technical variation you see from handling the samples. If the technical replicates are treated as biological replicates and differential analysis is performed, it will be seriously misleading.

I think both edgeR and DESeq2 doesn't estimate dispersion with 1 v 1 (which is right).

I would suggest maybe the following:

Run a fisher.test at every gene (methylated / unmethylated), get the odds ratio and p-value, rank this, but explain that the p-value is only for ranking, and will underestimate the variation.
A bit more refined, look for other studies that might be close to yours, get an estimate of the dispersion in their study and use a negative binomial distribution to again rank genes that show the biggest difference


Normalization using parallel R script I need to implement a parallel version of an R script because on one core it takes forever to execute. The function that takes the most of time is scran::quickCluster() applied to an integer sparse matrix (UMI_count) that is 26000x66000. The code is presented below:

library(scran)
library(scater)
sce <- newSCESet(countData=UMI_count)
clusters <- quickCluster(sce)
sce <- computeSumFactors(sce, cluster=clusters)

I did not find much on parallel versions of quickCluster() function, so will it actually execute in parallel if I set the number of cores to more than one? I am trying to do normalization of heterogeneous dataset with multiple cell types, that is why I am using quickCluster().

And I just need to speed it up to be able to compute it.

The scran::quickCluster method has a BPPARAM (i.e., BiocParallelParam), so if one provides something like

clusters <- quickCluster(sce, BPPARAM=BiocParallel::MulticoreParam())

it will use all available cores (a specific number of cores can also be given). See the intro doc to BiocParallel for info on other options.


If your data is in a dataframe you can make use of plyr. Using mapvalues() from this library, you can rename provide a list of keys and a list of values. Instead of explicitly listing the keys, we can just use unique(mat$Tumor_Sample_Barcode) [credit to @dc37]. I use the paste0 function here and the colon operator to generate the vector of A1,A2,...,A21.

Assign the vector produced by mapvalues() to the column and you're off!

install.packages("plyr")
mat$Tumor_Sample_Barcode <-
plyr::mapvalues(x=mat$Tumor_Sample_Barcode,
                from=unique(mat$Tumor_Sample_Barcode),
                to=paste0("A", 1:21))

I'd recommend keeping this column and adding a new column with your simplified sample names, but this is not what you requested.


Compute copy number from cases and controls I have some data on Copy Number Variation (SNP chip) for a population of samples.

In particular, I have a set of samples (considered as cases) which display a specific disease phenotype, and another set (considered as controls) which do not. The cases have not a matched-control. All the controls are taken as a random sample from the population that does not present the disease.

For both cases and controls I have the copy number of some regions. How can I compute the copy number change in cases, such as in GDC?

Should I use (as described in the previously linked page):

  performing tangent normalization, which subtracts variation that is found in a set of normal samples

and do you know of any tool to perform this computation?

The data I have is formatted in this way:

Sample_ID   chrom      start           end  CN
Sample11       19   11991477      12133823   1
Sample11        2   52260564      52431658   1
Sample12        7    5721757       5896192   3
Sample13       10    2269963       2473585   3

GISTIC does exactly this, and is very possibly the tool used in your link above.

The input files need to be in a different format than what you currently have, but if you're using a SNP6 array, there are guides out there that will tell you how to get the proper files from your .cel files.


R studio failed to install the latest development version of hdf5r from Github (win10) I install hdf5r following https://github.com/hhoeflin/hdf5r， it shows two methods:

When I use the second method to install from github,it showed many errors：

The downloaded source packages are in
‘C:\Users\dell\AppData\Local\Temp\RtmpI7TP9f\downloaded_packages’
√ checking for file 'C:\Users\dell\AppData\Local\Temp\RtmpI7TP9f\remotes8c87faaf3d\hhoeflin-hdf5r-cc278c1/DESCRIPTION'

* preparing 'hdf5r': (4.9s)
  √  checking DESCRIPTION meta-information ...
* cleaning src
* checking for LF line-endings in source and make files and shell scripts (676ms)
* checking for empty or unneeded directories (1.5s)
* building 'hdf5r_1.0.1.tar.gz'
  Warning: file 'hdf5r/cleanup' did not have execute permissions: corrected
  Warning: file 'hdf5r/configure' did not have execute permissions: corrected

* installing _source_ package 'hdf5r' ...
  checking hdf5 headers and libraries
  Environment variable LIB_HDF5 is not set. Attempting to download HDF5 files from repository.
  试开URL’https://github.com/mannau/h5-libwin/archive/master.zip'
  Content type 'application/zip' length unknown
  ** libs
  错误: (由警告转换成)this package has a non-empty 'configure.win' file,
  so building only the main architecture
* removing 'D:/soft ware/R-3.5.2/library/hdf5r'
  In R CMD INSTALL
  Error in i.p(...) :
  (converted from warning) installation of package ‘C:/Users/dell/AppData/Local/Temp/RtmpI7TP9f/file8c863b6795/hdf5r_1.0.1.tar.gz’ had non-zero exit status
  In addition: Warning messages:
  1: In untar2(tarfile, files, list, exdir) :
  skipping pax global extended headers
  2: In untar2(tarfile, files, list, exdir) :
  skipping pax global extended headers

But when I use the first command:install.packages("hdf5r"), it worked:

Can anyone help me? I search many times, but I could't find the answer.

Try running this:

source("https://bioconductor.org/biocLite.R")
biocLite("hdf5r")

Just a different method of installing packages besides github. It is more centred around bioinformatics packages as well so it is pretty useful.


Input normalization in ChIP-seq If I subtract input counts from ChIP counts (for every gene, since I have one peak per gene) I get negative values for most genes. This makes sense to me, because (as can be seen in the figure) input and IP have same sequencing depth but IP is very biased towards few genes with lots of counts.

Elsevier was informed that the business model which they offered does not meet the demands. How should I implement input normalization?

P.S. The subtraction idea is from Bioconductor

I strongly suggest that you not try to come up with your own package for this when things like CSAW already exist in bioconductor and provide a number of useful normalization options.

For visualization purposes, I think it's best to simply take the log2 ratio (e.g., using bamCompare from deepTools) or alternatively normalize by total coverage (possibly excluding blacklist regions, if applicable for your organism).


How can I get unmatched reads for defective genomes analysis using bwa and samtools? I am trying to follow the Beauclair et al. paper (free version here https://rnajournal.cshlp.org/content/24/10/1285.long) for identifying defective genomes using their DI-tector program.
According to figure 2
Before using the program itself, it is highly recommended to filter my raw read using bwa and then samtools
Unfortunately I have no idea where to start - I have browsed official bwa and samtools sites, I have installed them using the following commands:
sudo apt-get update -y
sudo apt-get install -y bwa

I have indexed my reference genome with the following command (as I read it is required before the analysis)
bwa index sequence.fasta

Then, I have installed samtools
cd ..
wget https://github.com/samtools/samtools/releases/download/1.9/samtools-1.9.tar.bz2
tar -vxjf samtools-1.9.tar.bz2
cd samtools-1.9
make
sudo apt install samtools

I wonder what will be the further steps to get these "UNMATCHED READS at the end of step i", since Di-tector program takes virus_reference genome (I have it) and these unmatched reads.
First of all I want to align my .fastq files vs GCF_000001405.39_GRCh38.p13_genomic.fna - what will be the output file? If I got it right - I need this output file to filter with samtools, am I correct?
If I am not mistaken, the output after samtools will be the .sam file that I will have to put into a DI-tector - correct?

If you want to keep track of the unmapped reads, I recommend you to use Bowtie2 instead of bwa.
Start by building an index:
bowtie2-build genome.fa genome.fa

Then map your reads, using arguments to save unmapped reads in a specific file like this:
bowtie2 -x genome.fa -1 reads_1.fq -2 reads_2.fq -S bowtie2_mapping.sam --un-conc-gz unmappedreads_%.fq.gz

It will write the alignment as well as unmappedreads_1.fq.gz and unmappedreads_2.fq.gz that you may use in DI-tector.


Mix globbing and wildcards when specifying rule input Consider the following scenario.

├── sample-alice
│   ├── sequence_1.fastq
│   ├── sequence_2.fastq
│   ├── ...
│   └── sequence_n.fastq
└── sample-bob
    ├── sequence_1.fastq
    ├── sequence_2.fastq
    ├── ...
    └── sequence_m.fastq

I'm trying to write a single rule in my Snakefile that will preprocess each sample generically. The rule will execute a single command that takes all of Fastq files associated with the sample simultaneously. The samples do not have the same number of Fastq files.

If it was a single sample, I could do something like this.

rule preprocess:
    input: glob('sample-alice/*.fastq')
    output: 'alice-clean.fastq'
    shell: 'mypreproccmd {input} > {output}'

On the other hand, if there was a single Fastq file for each sample, I could use a wildcard to write a generic rule like this.

rule preprocess:
    input: 'sample-{samp}/sequence_1.fastq'
    output: '{samp}-clean.fastq'
    shell: 'mypreproccmd {input} > {output}'

Is it possible to combine these two approaches to write a rule for handling multi-file samples generically?

A nice solution comes from a colleague of mine, as seen in this Snakemake workflow. The trick is to access the wildcards programmatically using an anonymous (lambda) function. In my example above, it would be implemented as follows.

rule preprocess:
    input: lambda wildcards: glob('sample-{samp}/*.fastq'.format(samp=wildcards.samp))
    output: '{samp}-clean.fastq'
    shell: 'mypreproccmd {input} > {output}'


Make ipyrad use cuda-enabled NVidia card on Ubuntu I want to use ipyrad on a new Ubuntu machine that has an NVidia Quadro K2000 card with 384 cores. One can configure ipyrad to run on a linux cluster.   Do I have any options to get ipyrad to access these cores?

I'm assuming that I will be able to install cuda drivers on the machine. If ipyrad can't be made to do this, is there another way to do de novo assembly of GBS data that will take advantage of this resource?

You would have to get CUDA Python set up on your system and likely modify ipyrad to actually use it (for example, by incorporating numba). So while this is possible, it will require significant effort on your part. In general, there are very few bioinformatics tools that support GPUs.


Resolution parameter in Seurat's FindClusters function for larger cell numbers In Seurats' documentation for FindClusters() function it is written that for around 3000 cells the resolution parameter should be from 0.6 and up to 1.2. I am wondering then what should I use if I have 60 000 cells? How to determine that?

Assuming you have an informative selection of variable genes from which you have constructed a number of useful PCs, I'd run a number of iterations with FindClusters() as described in the other answer, then choose a level which overclusters the dataset (for example, clusters that are visibly separate on a t-SNE or other dimensionality reduction plot should definitely have their own number):

seuratobject <- SetAllIdent(seuratobject, id='chosen.resolution')

Then run:

Seurat::BuildClusterTree()
Seurat::FindAllMarkersNode()

Assessing the cluster markers for each node will hopefully give you a good idea on which clusters should be combined. Then you can "combine" the clusters and re-label the cells using something like:

library(plyr)
cell.labels <- seuratobject@ident
cell.labels <- mapvalues(cell.labels,
  from=0:16, # cluster numbers
  to=c('A', 'B', 'C', 'C', 'D', 'E', 'E', ... )) # etc

seuratobject <- AddMetaData(seuratobject, cell.labels, 'Combined.clusters')

The usefulness of the clustering will very much depend on the selection of variable genes, therefore, depending on the (diversity of the) dataset, you will want to experiment with selection parameters or subset the dataset and repeat the above procedure.


Using large databases with BLAT I'm a computer scientist working with biologists at a small school that doesn't have dedicated bioinformatics staff. I apologize if I use incorrect terminology since I have limited bioinformatics background.

Our biology staff is studying a soil sample that has been sequenced. The resulting sequences (in the hundreds of millions) were analyzed using mg-rast, but they didn't seem to have a high confidence in the results.

I've been asked to help analyze these sequences. I setup a local instance of NCBI blast, but it was far too slow to analyze that many sequences. We moved on to using BLAT. They've requested we use the nt database from NCBI. I converted the nt FASTA file to 2bit file, splitting the files into small chunks.

gfServer seems to have limitations however. When I run gfServer, I get the following error while gfServer attempts to start.

  Exceeding 4 billion bases, sorry gfServer can't handle that.

Is it possible to use the full nt database using BLAT? Is there an alternative to BLAT that can handle that many sequences?

I'd recommend having a go with either OneCodex (commercial / Freemium) or Centrifuge.

Chris Mason's lab has recently looked at metagenomic classification and OneCodex gave the most accurate result.

Centrifuge uses Burrows-Wheeler lookups, but hasn't yet been compared by the Mason Lab.


how to obtain only one cluster's data (.csv file) in Seurat I was using Seurat to analysis single-cell RNA Seq. And I was interested in only one cluster by using the Seurat. Does anyone know how to achieve the cluster's data(.csv file) by using Seurat or any other method that can do this function? Thank you!

You can access data within the Seurat object using GetAssayData, and extract a list of cell names for the cluster you're interested in using WhichCells:

library(Seurat)

cells.use <- WhichCells(object = pbmc_small, ident = 1)
expr <- GetAssayData(object = pbmc_small, assay.type = "RNA", slot = "data")[, cells.use]
expr <- as(Class = 'matrix', object = expr)
write.csv(x = expr, file = "expression_cluster1.csv", quote = FALSE)


wget for links inside html pages I am trying to download a file from the following repository: https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR7276474

As you can see, there are several layers to the webpage. For example, clicking on the download tab doesn't change the URL and, the link is not a 'download link' per ce - where simply clicking the link automates a download. I have tried some of the answers on the forum, where they have advised using quotes and the operation:
wget "url/?target=link". This however does not work in the following instance.

Basically, I want to download the file labeled 'P1TLH.bam', in the download tab of the link provided:

If you want to download a file, you need to use the link to that file. Your original attempt, wget  https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR7276474 wouldn't work since that's a link to the trace page of the relevant run. If you want to download something else, just right click on the link (the one in the screenshot in your question), copy the URL and use that:

wget https://sra-download.ncbi.nlm.nih.gov/traces/sra64/SRZ/007276/SRR7276474/P1TLH.bam


Books on bioinformatics algorithms I'm looking for a book about bioinformatics algorithms, such as alignment, BLAST search, and variant calling.

I'm hoping reading about this subject will give me a deeper understanding of the foundations of bioinformatics, and I'm also interested in developing my own tools and algorithms for sequence and variant analysis. I studied applied mathematics and computer science at university, so I don't shy away from math-heavy books.

I've listed some candidates below.

Biological Sequence Analysis: Probabilistic Models of Proteins and Nucleid Acids, 1st edition, Durbin et al. 1998
An Introduction to Bioinformatics Algorithms, 1st edition, Jones et al. 2004
Algorithms in Bioinformatics: A Practical Introduction, Sung 2010

Durbin et al. is quite old, Sung is the newest.

I am looking for books with:

material (algorithms) on alignment, assembly, and variant calling if possible.
New books are preferred

Here are a couple books I'd recommend:

Dan Gusfield's Algorithms on Strings, Trees, and Sequences is a deep and wide treatment of aligning, searching, and processing strings, trees, and sequences.

Warren Ewens' Statistical Methods in Bioinformatics devotes a chapter to BLAST and the math underneath it.

Edit - Another book that may be useful for some mathematical underpinnings would be Michael S. Waterman's Introduction to Computational Biology.


how to have a blast xml file in a readable and understandable way I got a blast xml file after a DNA sequence search. The XML file is not very readable. What is the way to make it more understandable? This blast contains about twenty matches. Here is the xml file

The most human-readable blast output formats are 0-4, e.g. -outfmt 0 - which if I am not mistaken is the format used by web blast. However, these formats are a pain to automatically parse.
As Cowboy_Patrick pointed out, XML is more of a machine-readable format. The tabular output format ( -outfmt 6) is very commonly used, because it is easily parsed by unix command-line tools, and its structure also allows a somewhat easy human reading and understanding.


Combine Fastq by writeFastq is not working properly Currently I have 4 fastq file for each patient as WES data. For WES pipeline, I want to combine those fastq files. I am using writeFastq from ShortRead package in Bioconductor.But I think this function is not working properly.

The reason of my doubt is the the size the next step is quality check using FastQC tool. This tool is working on individual file but not on combined file and shows the error that ID line doesn't start with "@" . So I want to know what are the other effective tool for combining multiple fastq files without any loss of information.

Let me know if I am doing something wrong. Thanks

[EDIT-1] : The codes for combining multiple fastq files into one is as follows:

reference_fastq = file.path(paste(getwd(),"refs",sep = "/", collapse = " "), fsep = ".fasta")
fastqDir <- file.path(dataDir, fsep = ".fastq")
outfile_temp <- paste(dirname(dataDir),basename(fastqDir),    paste(basename(fastqDir),"combined","fastq", sep = "."),sep = "/", collapse = " ")
fls <- list.files(fastqDir, "fastq$", full=TRUE)
for (fl in fls) {
  fq <- readFastq(fl)
  writeFastq(fq,outfile_temp, mode="a")
}

First four lines of my fastq file (one of 4 indivual file) are as follows :

@D00239:185:C9KELANXX:2:2201:1088:2107 1:N:0:TCCTGAGC
GTACAAGAAGCATAACATCAGCATGTGTTTTGGTGAGAACCTTAGGAAGTTTGCAATCATGACAGAAGGTGAAGG
+
/<BBBFFFBFFFFFFFFFFFBF<FFBF<B<FFFFFFFFFFB/B/<////<<<FF<B<FF<FFFFFFFFF/FFFFF

First four lines of combines fastq file are as follows:

@D00239:185:C9KELANXX:2:2201:1088:2107 1:N:0:TCCTGAGC
GTACAAGAAGCATAACATCAGCATGTGTTTTGGTGAGAACCTTAGGAAGTTTGCAATCATGACAGAAGGTGAAGG
+
/<BBBFFFBFFFFFFFFFFFBF<FFBF<B<FFFFFFFFFFB/B/<////<<<FF<B<FF<FFFFFFFFF/FFFFF

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

Thanks.

You can drastically simplify your effort because FASTQ files are simple text files; therefore, standard UNIX test file tools work:

cat refs/*.fastq > combined.fastq

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

Apart from this, your R code contains several errors. For instance, your variable names change throughout, and you’re using the fsep parameter wrongly. In fact it’s unclear to me what the code is supposed to do. At an approximation (and untested), the following should work:

fastq_dir = file.path(getwd(), "refs")
fastq_files = dir(fastq_dir, pattern = "\\.fastq$", full.names = TRUE)
outfile = file.path(dirname(fastq_dir), "combined.fastq")

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

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


Trim Seurat object to contain expression info only for selected genes I'm working with some large Seurat objects (MOCA, MCA, Tabula Muris) studying gene coexpression, and I'm running into memory issues.

Is it possible to remove all genes in a Seurat object that are not in a specified gene array? Ideally I'd like to be working with a small R matrix with say, 250k rows (cells) and 10-50 columns (gene expression data) which will enable fast processing and comprehensive analysis using other packages.

A different approach if you are using Seurat3, is DietSeurat(). It allows you to diet the object by removing the components that you don't need. For example you can keep the normalised/scaled matrix and remove the raw counts. This approach could reduce space and memory usage, while keeping all your genes in place.


Plot a circos plot to show the consistency between 2 samples I have called SNP and INDEL in two matched samples by strelka and extract this information from .vcf file and I have these

CHROMOSOME  POS REF ALT SAMPLE
1   928006  G   A   t_005
1   1649842 G   T   t_005
1   2020408 G   A   t_005
1   2031677 T   A   t_005

and

CHROMODOME  POSITION    REF ALT SAMPLE
1   14115   A   T   o_005
1   541052  T   C   o_005
1   1088123 T   G   o_005
1   1232501 A   G   o_005

I need to show if and how extend the mutations between these two samples are consistent by a circos plot but really I don't know how people do that

Something like this

or

A circos plot is most likely not the most appropriate solution here. What I would suggest is a confusion matrix, of which you can find an example here:

For every variant in your vcf you'll add a number in this matrix. One sample is the columns, the other is the lines. If your variant is homozygous in both, then you add in that square +1 (the cell with 5845 in the example).

A perfect concordant sample pair will have only variants on the diagonal.

Here is some python code to get such a matrix. It uses cyvcf2 and pandas, and expects as input a vcf file with both samples.

from argparse import ArgumentParser
from cyvcf2 import VCF
import pandas as pd

def main():
    args = get_args()
    confusion_matrix(args.vcf)

def confusion_matrix(vcff):
    """
    First level of the dict is the "first" call, second level is the "second" sample
    0: hom_ref
    1: heterozygous
    2: unknown/nocall
    3: hom_alt
    """
    zygosities = {0: {0: 0, 1: 0, 2: 0, 3: 0},
                  1: {0: 0, 1: 0, 2: 0, 3: 0},
                  2: {0: 0, 1: 0, 2: 0, 3: 0},
                  3: {0: 0, 1: 0, 2: 0, 3: 0},
                  }
    for v in VCF(vcff):
        zygosities[v.gt_types[0]][v.gt_types[1]] += 1
    zygs = [2, 0, 1, 3]
    df = pd.DataFrame(index=zygs, columns=zygs)
    for tr in zygs:
        for te in zygs:
            df.loc[tr, te] = zygosities[tr][te]
    df.columns = ['nocall', 'hom_ref', 'het', 'hom_alt']
    df.index = ['nocall', 'hom_ref', 'het', 'hom_alt']
    print(df)

def get_args():
    parser = ArgumentParser(description="Create confusion matrix of SNV calls")
    parser.add_argument("vcf", help="vcf containing two samples")
    return parser.parse_args()

if __name__ == '__main__':
    main()