Extracting expression data in r from GEO and got a S4 object I want to open the expression data GSE9838 dataset from GEO using GEOquery R package.

gset <- getGEO(filename="GSE9838_family.soft.gz")

The result of this code said it's S4 object of class GSE.

How can I get the expression data of GSE9838?

Well, if you use a SOFT formatted files:

#You can read your file
gds <- getGEO(filename="GSE9838_family.soft.gz")

# and check it
class(gds)
show(gds)

#then you usually would try to use  an exprs() accessor to retrieve the expression matrix:
exprs (gds[[1]])

# but in this case you will get an error revealing that it is a GSE datastructure
# eg. is composite structure having GSM and GPL object lists:
names(GSMList(gds))
names(GPLList(gds))

# now GPLList is probably not interesting in this case
# but you can look at the GSM-objects from GSMList:
first <- GSMList(gds)[[1]]
first
class(first)
Meta(first)

# at this point if you read the results you´ll notice that
#    title
#    [1] "X_intercept"

# and
#   An object of class "GEODataTable"
#   ****** Column Descriptions ******
#   [1] Column      Description
#   <0 rows> (or 0-length row.names)
#   ****** Data Table ******
#   data frame with 0 columns and 0 rows

So here you would probably wonder if the data you want is these objects. If you go back go GEO page for GSE9838, you can find that "Raw data provided as supplementary file. Processed data is available on Series record". So one may want to check series record instead:

# you can query the GSE9838 like this:
data = getGEO("GSE9838")

# then you can check what you´ve got
names(data)
#[1] "GSE9838_series_matrix.txt.gz"
class(data$GSE9838_series_matrix.txt.gz)

# and check what´s inside:
exprs(data$GSE9838_series_matrix.txt.gz)
# but apparently it just shows one header row.

So it wasn´t there either. And if you check the size of the file, you´ll find that it is just 4.8 kb. My educated guess would be that the expression values are not in the dataset itself but probably in some of the supplementary files?

If you want to download them, you can:

filePaths = getGEOSuppFiles("GSE9838")
filePaths
untar("GSE9838/GSE9838_RAW.tar", exdir = "GSE9838/CEL")
list.files("GSE9838/CEL")

...what will allow you to see all the CEL files with the data.


Are soft-clipped bases used for variant calling in samtools + bcftools? If there are soft clipped base pairs specified in the CIGAR string for a read in a SAM/BAM file, will these be used for variant calling in a samtools + bcftools workflow?

The GATK HaplotypeCaller, for example, has an explicit option --dontUseSoftClippedBases for whether to use soft clipped bases. The samtools documentation does not mention clipped bases.

No, samtools (and therefore bcftools) does not use soft-clipped bases. You can quickly confirm this by using either samtools depth or samtools mpileup to look at a region with a soft-clipped alignment. You'll note that the soft-clipped region isn't used in the depth/pileup (both tools use the same underlying code, so it doesn't matter which you use). If you're curious, samtools ignores soft-clipped bases because it's based on making a per-base stack of alignments covering each position. In the BAM format, alignments are sorted and assigned to bins according to their start/end positions, which won't include soft-clipping. Consequently, when samtools is making the pileup it won't even see the alignments that would overlap a given base if soft-clipped bases were included.

This then sort of begs the question of what GATK's HaplotypeCaller is doing differently. There, regions in the genome are essentially assembled in a small de Bruijn graph, which allows for soft-clipped bases around indels to then be resolved, given that the graph would start/end a little-way on past each side of indels. This is also why you don't need to do indel realignment with the HaplotypeCaller (this was needed in the old UnifiedGenotyper).

Edit: For more details regarding the HaplotypeCaller, see this nice page on GATK's website, which goes into much more detail than I did here.


Plink 1.9 equivalent of --score sum in Plink 2.0 I have previously performing Allelic scoring in Plink 1.9 using:

plink \
--score sum sumstats.txt\

However, it seems this has changed in Plink 2.0. What is the equivalent function for Allelic scoring?

plink2 --score cols=+scoresums

See https://www.cog-genomics.org/plink/2.0/general_usage#colset and https://www.cog-genomics.org/plink/2.0/score for details.


How many reads has my sequencing run produced on minion? I am running a 48 h sequencing protocol on a FLO-MIN107, on MinION, with DNA library-prepped with SQK-RAD004. MinKNOW is installed on my Mac and controlling the MinION. I have found the following histogram in MinKNOW. What do the numbers mean? In particular, how many reads have been produced? I understood that the events are consecutive electrical signals and are then transformed into reads? Do I have 4029 reads? The run has started 24 hours ago. 4029 reads surely is really low, isn't it? What does 4,120,722 mean?

The Platform QC of the flow cell was as follows:

group 1: 480
group 2: 374
group 3: 201
group 4: 59
total = 1114

That's indeed the number of reads and that's quite low. How was your pore occupancy (number of pores sequencing) and your flow cell QC (number of pores good enough for sequencing)? How was your library concentration?

Events is an approximation of nucleotides: it's a guess based on the current signal. As far as I know if you multiply this by 1.8 you get approximately your throughput in nucleotides.

The 4,120,722 is the number of events in the currently selected (mouse hovering?) bin of the histogram.


How to create a nexus file for BEAST? I am learning to use BEAUTI and BEAST software. I followed their tutorial and was able to run example data without any problem. Now, I want to do phylodynamic analysis on my data. I have downloaded 70 picornavirus sequences from the gene bank. My questions are :

1) What tool do I have to use to create Nexus file that is used by BEAUTI ?

2) How do I enter the information about the date (year ) the virus was isolated and location information which is required for rates and dates (evolutionary) analysis in BEAST, spatial analysis by SpreaD3 ?

Its easy,

Download the data from NCBI in fasta format (its an option in one of the pull down headers).
Load it into the latest version of MEGA, version 7(?), perform the alignment therein.
Export from MEGA as Nexus format and upload it into Beast

You could of course perform the alignment in Clustalx, Clustalo or Muscle or T-Coffee, (hmmm I could write an alignment program just to name it T-4two (sorry British joke)). If the output option of "nexus" was not present in the alignment program of your choice, then again upload the alignment into MEGA and export it as nexus.

Anyway, I could go into specific detail about the various points and clicks therein, but I think its easy to figure out. Most phylogeny programs will do an automatic fasta->nexus conversion, but anyway.

Its easy, Beast is not trivial to use properly however and ensuring convergence is not at all easy, in addition the multiple simulation a priori models, whilst powerful, give different answers. To try and explain this a different way, it is easy to get an answer in Beast, but difficult to get a robust answer from it and thats the bit I would give serious consideration about. It does amaze how widespread the program is (I guess its the nice GUI), given how complicated it is to obtain robust result.

Just to add you could try the online tool here for the fasta to nexus conversion:  . I've not tried it, its just an after thought. If you do Python, Biopython will do it, but its just as easy to code.


Finding discriminating positions for a sequence alignment column I don't have any theoretical background in biology, so forgive me if my question is a bit.. well..dumb. I'm trying to use a Monte Carlo approach to find the discriminating position of a given sequence alignment column. So essentially, given a sequence alignment, I want to find the position at which we can distinguish active and inactive conformations.

So, here is my code for the Monte Carlo method:

def MonteCarlo(S, max_iter, threshold):
    """
    Computes the discriminating position for a sequence column

    Input: S = sequence column, max_iter = max iterations, threshold = threshold of acceptance for boltzmann distribution (between 0 and 1)
    Output: 2 classes
    """
    n = len(S)
    S_copy = S.copy()
    num_class1 = int(0.30*n) #number in class1
    print(num_class1)
    for _ in range(num_class1): #intialisation
        rand_position = random.randint(0,len(S)-1)
        if S_copy[rand_position] == 'X':
            rand_position = random.randint(0,len(S)) #reselect if used position
        S_copy[rand_position] = 'X'

    #Monte Carlo with Insertions and Deletions

    s_entropy = entropy(S)
    n_iter = 0 #counter for iterations
    class_1 = dict() #class 1 determined by X's
    class_2 = dict() #class 2 is rest of S
    for i in range(len(S)):
        if S_copy[i] == 'X':
            class_1[i] =S[i]
        else:
            class_2[i] =S[i]

    res = score(S, class_1, class_2)

    while n_iter < max_iter:
        r = random.random()

        if r < 0.25: #25% chance to do an insertion
          (c1,c2,pos) = Insertion(S, class_1, class_2)

        elif r < 0.5: #25% chance to do a deletion
          (c1,c2,pos) = Deletion(S, class_1, class_2)

        else: # 50% chance to do a swap
          (c1,c2,pos) = Swap(S, class_1, class_2)

        res2 = score(S, c1, c2)

        #Boltzmann Distribution coefficient
        delta_S = abs(res2 - res)
        lambda_, N = 0.2, 19
        size = 6000
        ind1 = random.randint(0,size - 1)
        rv = boltzmann.rvs(lambda_, N, size=size)/N
        beta = rv[ind1] #pick float between 0 and 1 with boltzmann distribution
        if  (res2 > res) or (beta*delta_S > threshold): #if better score i.e. higher entropy gain
            res = res2
            class_1 = c1
            class_2 = c2

        n_iter += 1

    return class_1, class_2, pos

Here I'm using the Boltzmann distribution as a probability measure. All the other functions (swap, insertion, deletion, entropy) calculate pretty much what the name says (if you need more clarification on my side, I'd be happy to provide it!)

Let me explain what my code is trying to do: Given a sequence alignment (as a list) say ['M','M','M',...,'L','L','L'] I first set a bunch of random positions and mark them out as 'X' (the number of Xs are determined by num_class1). Next, in the while loop I do swaps,insertions, and deletions with the goal of obtaining a better score. This meansI'm adding, removing, or swapping positions between class1 and class2. Eventually, my code returns two classes and the discriminating position.

While my code works, I don't think it's giving the right result. When I use a base case example of a list with two letters 'M' and 'L with a clear discriminating position, I get random results.

Could someone explain theoretically (or via code) how one can find the discriminating position using a Monte Carlo method as above?

Thanks!

Without asking further questions, my suspicion would be how indels (gaps) are counted in your code against the aligment. Its the classic place for alignment mis-position errors. In other words, based on my experience your calculation probably performs fine, but linking the calculated alignment position against the actual observed alignment position is classic area of logical errors and needs careful monitoring of variables using code breakpoints within your loops.


10X cellranger count, [error] The chemistry was unable to be automatically determined While running cellranger count using the following .slurm file:

#SBATCH -p standard
#SBATCH -A overall
#SBATCH --time=48:00:00
#SBATCH --mail-type=end
#SBATCH --mail-user=nv4e@virginia.edu
#SBATCH --output=generate_count_%j.out
#SBATCH --error=generate_count_%j.err
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=20
#SBATCH --mem=120000

# Run driver program
FASTQS="/scratch/nv4e/kipnis/AHN3LLBGX5/outs/fastq_path/HN3LLBGX5"
TR="/scratch/nv4e/kipnis/refdata-cellranger-mm10-1.2.0"

cellranger count --id=17R_1 --transcriptome=${TR} --fastqs=${FASTQS} --sample=17R --expect-cells=2000

I got the error:

>2018-04-28 20:04:37 [runtime] (progress)
ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMI
STRY.fork0.chnk0: Indexing genome...
2018-04-28 20:04:49 [runtime] (progress)        ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMISTRY.fork0.chnk0: Building transcriptome...
2018-04-28 20:04:58 [runtime] (progress)        ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMISTRY.fork0.chnk0: Building kmer index...
2018-04-28 20:05:25 [runtime] (progress)        ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMISTRY.fork0.chnk0: Mapping reads...
2018-04-28 20:05:28 [runtime] (chunks_complete) ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMISTRY
2018-04-28 20:05:28 [runtime] (run:local)       ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMISTRY.fork0.join
2018-04-28 20:05:31 [runtime] (join_complete)   ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMISTRY

2018-04-28 20:05:34 [runtime] (ready)           ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.SETUP_CHUNKS
2018-04-28 20:05:37 [runtime] (split_complete)  ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.SETUP_CHUNKS
2018-04-28 20:05:37 [runtime] (run:local)       ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.SETUP_CHUNKS.fork0.chnk0.main
2018-04-28 20:05:40 [runtime] (failed)
ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.SETUP_CHUNKS

>[error] The chemistry was unable to be automatically determined. This can happen if not enough reads originate from the given reference. Please verify your choice of reference or explicitly specify the chemistry via the --chemistry argument

What is going wrong here? How could I correct it?

The issue was that I supplied wrong indexes during fastq files generation, however, they were generated successfully, and only the next step failed with not that informative error.


Using Beagle 4.1 for phasing and IBD - chromosomes I am trying to use beagle 4.1 for phasing and IBD analyses. I was wondering if I need to run the analysis by chromosome or I can put them together in the same vcf file.

OK, I ended up contacting one of the authors of the program, who kindly answered to my question:

  Either approach should work,

[i.e. running the analysis by chromosome vs joining the data together]

  but I suggest dividing by chromosome because that approach has had
  more testing.

I hope this might help other users, I saw that several questions regarding Beagle went unanswered on biostars and seqanswers. I guess it is a pretty specialized topic.

Best,

Max


What process and input data is required for a cellranger reference transcriptome? I'm analysing single-cell RNA-Seq data using the 10X Genomics cellranger platform. While they provide reference data for Mouse and Humans, other species require a reference "transcriptome" to be generated with the cellranger mkref tool. However, I've found it to be very temperamental giving errors such as:

gffread genome_GFF3_genes.gff -T -o genome_genes.gtf

cellranger mkgtf genome_genes.gtf genome_genes.filtered.gtf \
                   --attribute=gene_biotype:protein_coding \
                   --attribute=gene_biotype:lincRNA \
                   --attribute=gene_biotype:antisense \
                   --attribute=gene_biotype:IG_LV_gene \
                   --attribute=gene_biotype:IG_V_gene \
                   --attribute=gene_biotype:IG_V_pseudogene \
                   --attribute=gene_biotype:IG_D_gene \
                   --attribute=gene_biotype:IG_J_gene \
                   --attribute=gene_biotype:IG_J_pseudogene \
                   --attribute=gene_biotype:IG_C_gene \
                   --attribute=gene_biotype:IG_C_pseudogene \
                   --attribute=gene_biotype:TR_V_gene \
                   --attribute=gene_biotype:TR_V_pseudogene \
                   --attribute=gene_biotype:TR_D_gene \
                   --attribute=gene_biotype:TR_J_gene \
                   --attribute=gene_biotype:TR_J_pseudogene \
                   --attribute=gene_biotype:TR_C_gene

cellranger mkref --genome="reference-name" \
                 --fasta="transcriptome.fasta" \
                 --genes="genome_genes.filtered.gtf"

  Fatal INPUT FILE error, no valid exon lines in the GTF file: reference_name/genes/genes.gtf
  Solution: check the formatting of the GTF file. Most likely cause is the difference in chromosome naming between GTF and FASTA file.

  subprocess.CalledProcessError: Command '['STAR', '--runMode', 'genomeGenerate', '--genomeDir', 'reference_name/star', '--runThreadN', '1', '--genomeFastaFiles', 'reference_name/fasta/genome.fa', '--sjdbGTFfile', 'reference_name/genes/genes.gtf', '--limitGenomeGenerateRAM', '17179869184', '--genomeSAsparseD', '1', '--genomeSAindexNbases', '12', '--genomeChrBinNbits', '10']' returned non-zero exit status 104

  subprocess.CalledProcessError: Command '['STAR', '--runMode', 'genomeGenerate', '--genomeDir', 'reference_name/star', '--runThreadN', '1', '--genomeFastaFiles', 'reference_name/fasta/genome.fa', '--sjdbGTFfile', 'reference_name/genes/genes.gtf', '--limitGenomeGenerateRAM', '17179869184', '--genomeSAsparseD', '1', '--genomeSAindexNbases', '12', '--genomeChrBinNbits', '10']' returned non-zero exit status -11

What input data is required for the cellranger to construct a reference? What format do the fasta and gtf files need be in? Is it possible to use a "transcriptome" rather than a "genome" as a reference here? For example, I have a transcriptome reference with fasta names >name for each transcript rather than each chromosome.

Your problem is caused by using the transcriptome fasta file rather than the genome fasta file. You've already given it transcriptome information with genome_genes.filtered.gtf, it needs the genomic sequence to go along with that.

As an aside, your GTF filtering command is unlikely to be useful unless you have a Gencode or Ensembl GTF for human. Have a look at the gene_biotypes in your file and decide for yourself what (if anything) needs to be filtered. For example, I made a cellranger reference last night for dm3 (the old fruit fly genome) with:

cellranger mkgtf annotation.gtf genes.gtf \
    --attribute=gene_biotype:protein_coding \
    --attribute=gene_biotype:ncRNA \
    --attribute=gene_biotype:snRNA \
    --attribute=gene_biotype:snoRNA \
    --attribute=gene_biotype:pre_miRNA

That should be customized for relevance. It's quite possible that your GTF doesn't even have gene_biotype attributes, that's specific to Ensembl I think.


Functional Annotation vs Functional Enrichment--which one should I use for Network Analysis? I am currently working on a network pharmacology project where I am supposed to find out the possible diseases our drug might have an effect on. I am not a bioinformatics student, so my knowledge in the subject matter is very limited. I have all the targets and I am supposed to get the pathways associated with the targets. I used STRING for the purpose.

I got both KEGG functional annotation and enrichment for my list of targets.
I have to build a two-mode network

Question I am not very clear about

the difference between enrichment and functional annotation and
which data should be used for the next step.

The answer to your main doubt about "annotation vs. enrichment" is that you need both, but you already have one annotation.
You can only perform a functional enrichment analysis after you have a functional annotation. In your question, you state that you have KEGG annotation for your targets. That is it, knowing which target "belongs" to each KEGG module is an example of functional annotation of your targets. Several tools can annotate genes, proteins, etc. Knowing which protein is related to each drug or disease is also a functional annotation.
Now, knowing which targets are significantly important can be achieved by functional enrichment analysis.
For the network, both information can be used together. For example, you could link all elements to their functional annotation and use the available enrichment score as a quantitative parameter of importance for calculating correlations.


Can I use DRIMseq to model interaction effect of variables? I am trying to use DRIMseq for DTU with 2 treatments on two different strains of mice. In my previous analysis with genes only, I had found a strong interaction effect  of strain and treatments. Is it possible to model that with DRIMseq?

Also when I try to add strain term and use it as the batch effect as described in the tutorial, I get the following error:

! Using a subset of 0.1 genes to estimate common precision !
Error in optimHess(par = par, fn = dm_lik_regG, gr = dm_score_regG, x = x,  :
  non-finite value supplied by optim
In addition: There were 50 or more warnings (use warnings() to see the first 50)

Please see the sessionInfo() below.

R version 3.4.0 (2017-04-21)
Platform: x86_64-redhat-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so

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

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

other attached packages:
 [1] DRIMSeq_1.6.0                            pheatmap_1.0.8                           ggplot2_2.2.1
 [4] edgeR_3.16.5                             limma_3.30.13                            TxDb.Mmusculus.UCSC.mm10.knownGene_3.4.0
 [7] GenomicFeatures_1.26.4                   AnnotationDbi_1.36.2                     Biobase_2.34.0
[10] GenomicRanges_1.26.4                     GenomeInfoDb_1.10.3                      IRanges_2.8.2
[13] S4Vectors_0.12.2                         BiocGenerics_0.20.0                      tximport_1.2.0

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.12               RColorBrewer_1.1-2         plyr_1.8.4                 pillar_1.2.1               compiler_3.4.0
 [6] XVector_0.14.1             bitops_1.0-6               tools_3.4.0                zlibbioc_1.20.0            biomaRt_2.30.0
[11] digest_0.6.12              bit_1.1-12                 gtable_0.2.0               RSQLite_2.0                memoise_1.1.0
[16] tibble_1.4.2               lattice_0.20-35            pkgconfig_2.0.1            rlang_0.2.0                Matrix_1.2-10
[21] DBI_0.7                    stringr_1.3.0              rtracklayer_1.34.2         Biostrings_2.42.1          locfit_1.5-9.1
[26] bit64_0.9-7                grid_3.4.0                 XML_3.98-1.10              BiocParallel_1.8.2         magrittr_1.5
[31] reshape2_1.4.3             blob_1.1.0                 MASS_7.3-47                scales_0.5.0               Rsamtools_1.26.2
[36] GenomicAlignments_1.10.1   SummarizedExperiment_1.4.0 colorspace_1.3-2           labeling_0.3               stringi_1.1.7
[41] lazyeval_0.2.1             munsell_0.4.3              RCurl_1.95-4.10

One way to see the interaction effect would be including it in the design matrix as explained in the limma tutorial. Something on the line of:

design_full <- model.matrix(~ treatment*strain, data = samples(d))

Or you could build it yourself in a matrix like

 (Intercept) strainA_treatmentA strainA_treatmentB  strainB_treatmentA strainB_treatmentB
       1                      0                  1                  0                   0
       1                      1                  0                  0                   0

About the error, it seems like the data is not in the right format, you have either NA or Inf somewhere in the matrix. Check it and rethink the processing you made.


How do I lift GWAS results to hg38? There are some questions about lifting between reference builds, e.g. this one. But there doesn't appear to be a question about lifting a GWAS results file to a new reference build (except off-site). GWAS results files are not VCFs or plink files, so they aren't covered by the usual lifting tools (such as CrossMap or LiftOver). How do I lift a GWAS results file to hg38?
For instance, the Lee et al 2018 GWAS on education has this file:
user@computer:/home/user/Lee_2018$ head GWAS_EA.to10K.txt
MarkerName  CHR POS A1  A2  EAF Beta    SE  Pval
rs9859556   3   49455986    G   T   0.6905  -0.02901    0.00151 4.61e-82
rs7623659   3   49414791    T   C   0.3095  0.02899 0.00151 6.05e-82
rs11917431  3   49644012    C   T   0.6973  -0.0292 0.00152 7.40e-82
rs1873625   3   49666964    C   A   0.6973  -0.02916    0.00152 8.10e-82
rs11921590  3   49644193    T   C   0.6973  -0.02916    0.00152 1.17e-81
rs2352974   3   49890613    T   C   0.5 -0.02569    0.00141 6.35e-74
rs7029718   9   23358495    G   A   0.5646  -0.0244 0.00143 7.82e-65
rs7868984   9   23357826    C   T   0.4388  0.02437 0.00143 8.21e-65
rs11793831  9   23362311    T   G   0.4354  0.02418 0.00143 9.34e-64

It turned out to be easier than expected. You figure out which dbSNP version they used. dbSNP is the classification system for snps based on the rsid format. Pretty much every GWAS publishes their results with the rsids for the variants, as done above as well (MarkerName). So you download one of the dbSNP files from the official site. In this case, I am using version 142 since this was used by the 1000 genomes dataset.
After this, it is a matter of joining the tables. This can be done easily in R (you need about 15GB free memory), we load the dbSNP database into memory. Since we don't need all the columns, we skip some of them to save memory:
#library
library(tidyverse)

#cols
vcf_cols = cols_only(
  `#CHROM` = col_character(),
  `POS` = col_number(),
  `ID` = col_character(),
  `REF` = col_character(),
  `ALT` = col_character()
  )

#this is from the 1kg
#ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/other_mapping_resources/
dbsnp142 = readr::read_tsv(
  "ALL_20141222.dbSNP142_human_GRCh38.snps.vcf"),
  comment = "##",
  col_types = vcf_cols
  )

GWAS_EA = read_tsv(GWAS_EA.to10K.txt,
                   col_types = cols(chrpos37 = col_character(),
                                    chrpos38 = col_character()))

#make chrpos
dbsnp142 %<>% mutate(
  chrpos38 = str_c((CHROM %>% str_replace("chr", "")), ":", POS)
)

#count overlap
(GWAS_EA$MarkerName %in% dbsnp142$ID) %>% table2()

#mutate and merge
GWAS_EA = GWAS_EA %>%
  mutate(
    chrpos37 = str_c(CHR, ":", POS)
  ) %>%
  left_join(dbsnp142 %>% select(ID, chrpos38), by = c("MarkerName" = "ID"))

#update file on disk
GWAS_EA %>% write_tsv("GWAS_EA.to10K.txt")

The two new columns, chrpos37 and chrpos38, have the chr:pos variant names in both versions.


Tools for modelling and visualizing growth of cells I am looking for a program/GUI app/package that would help me do simulations for cell growth. I have a microscopic mathematical model, and the scenario is basically the following: I start with a rod-shaped cell characterized my its center and direction angle, it grows until it reaches a certain length, it divides into two identical cells whose angles change randomly and the process continues.

Is there any program or package that I can use for such a situation? I have found several on GitHub, but they were mostly for very complicated biological processes. I use Python/Julia/Matlab and can code, but I was wondering whether there are any tools for such simpler scenarios.

If the scenario is simpler, and you can't find anything by searching, I would consider just writing your own.

The modelling can be done in any language but if compute time is an issue - for example you want to run a simulation with many cells or many simulations - then I would use a language that can support that such as Julia or NumPy in Python. In fact Julia is a nice language for this type of thing as you end up with readable code that runs fast.

For visualisation in Julia specifically you could use the animation functionality of Plots.jl, which is basic but there, or could try Makie.jl. See the Molly.jl documentation for some examples of visualising particle-type processes using Plots.jl.


Setting Contrast for DESeq2 results This question was also asked on BioStars
I am trying to recreate this heatmap. How would I compare the four variables: Ly49+/- and MOG/MOGSP to get a single heatmap? Thank you for helping me through this.
I need to compare both Ly49+ versus Ly49− and MOG versus MOG plus SP. Column data in at the bottom. The function I used for DESeq2 is,
dds = DESeqDataSetFromMatrix(countData=countData,

colData=colData,

design= ~ cell_treatment)

I am not sure if the design is right. For the results, how should I write the contrasts? I am very confused and would really appreciate any help!
DataFrame with 11 rows and 2 columns

cell treatment

<character> <character>

MOGSP_2 SPtetramer+CD8+TCell MOGSP

MOGSP_3 SPtetramer+CD8+TCell MOGSP

MOG_3 SPtetramer+CD8+TCell MOG

MOG_4 SPtetramer+CD8+TCell MOG

MOG_5 SPtetramer+CD8+TCell MOG

Ly49N_1 Ly49-CD8+TCell MOGSP

Ly49P_1 Ly49+CD8+TCell MOGSP

Ly49N_2 Ly49-CD8+TCell MOGSP

Ly49P_2 Ly49+CD8+TCell MOGSP

Ly49N_3 Ly49-CD8+TCell MOGSP

Ly49P_3 Ly49+CD8+TCell MOGSP

Thank you!!

Your table is not very clear. If I understood it correctly, there are three groups of t-cells, SPtetramer+, Ly49- and Ly49+ and within it there is a MOGSP vs MOG comparison. I would set colData like this:
                   name             Ly treatment
1  SPtetramer+CD8+TCell SPtetramerplus     MOGSP
2  SPtetramer+CD8+TCell SPtetramerplus     MOGSP
3  SPtetramer+CD8+TCell SPtetramerplus       MOG
4  SPtetramer+CD8+TCell SPtetramerplus       MOG
5  SPtetramer+CD8+TCell SPtetramerplus       MOG
6        Ly49-CD8+TCell      Ly49minus     MOGSP
7        Ly49+CD8+TCell       Ly49plus     MOGSP
8        Ly49-CD8+TCell      Ly49minus     MOGSP
9        Ly49+CD8+TCell       Ly49plus     MOGSP
10       Ly49-CD8+TCell      Ly49minus     MOGSP
11       Ly49+CD8+TCell       Ly49plus     MOGSP

I generate an example dataset:
counts = matrix(rnbinom(11000,mu=100,size=1),ncol=11)

And in the design you specify these two factors, meaning gene expression can be affected by these two conditions:
dds = DESeqDataSetFromMatrix(counts,colData=colData,~Ly +treatment)
dds = DESeq(dds)

You can compare them like this:
head(results(dds,c("Ly","Ly49plus","Ly49minus")))
log2 fold change (MLE): Ly Ly49plus vs Ly49minus
Wald test p-value: Ly Ly49plus vs Ly49minus
DataFrame with 6 rows and 6 columns
          baseMean     log2FoldChange            lfcSE               stat
         <numeric>          <numeric>        <numeric>          <numeric>
1 91.4458114779915  0.712987165797688 1.28291188260395  0.555756927241586
2 81.7945234843954 -0.418266530066077 1.14378290088875 -0.365686993345568
3 94.4351681904941 -0.775592496878288 1.32965355129901 -0.583304196886903
4 107.484335203298  -1.36694885973712 1.38595911743023 -0.986283680770933
5 61.2205696982841 -0.760219864564472 1.16779888688769 -0.650985262188887
6 95.5908168437137  0.497427543096335 1.12768919469662   0.44110340458672
             pvalue              padj
          <numeric>         <numeric>
1 0.578377034355333 0.977537706056193
2    0.714598652706 0.992321035598301
3 0.559688537757548 0.977537706056193
4 0.323993925580938 0.954772990635094
5 0.515056000552832 0.977537706056193
6 0.659138138867294 0.992321035598301

To compare groups within another factor:
results(dds,c("treatment","MOGSP","MOG"))
log2 fold change (MLE): treatment MOGSP vs MOG
Wald test p-value: treatment MOGSP vs MOG
DataFrame with 6 rows and 6 columns
          baseMean     log2FoldChange            lfcSE               stat
         <numeric>          <numeric>        <numeric>          <numeric>
1 91.4458114779915 -0.875905107403943 1.43007567177694 -0.612488642867122
2 81.7945234843954  0.595588650280833 1.28109080878515  0.464907441530728
3 94.4351681904941 0.0255578400583986 1.48592401497315 0.0171999643325371
4 107.484335203298  -1.46652057431159 1.55285829334496 -0.944400774105796
5 61.2205696982841  -1.27442596494316 1.30489822606005  -0.97664778715433
6 95.5908168437137  -2.35402796630456 1.26608624796739  -1.85929510733078
              pvalue              padj
           <numeric>         <numeric>
1  0.540214509839105 0.944881849859735
2  0.641997741235223 0.973651054534338
3  0.986277090644156 0.994944323499595
4  0.344964886056295  0.87421063696937
5  0.328743552427222                NA
6 0.0629853201285643 0.525889480109579

The colData i used:
structure(list(name = structure(c(3L, 3L, 3L, 3L, 3L, 1L, 2L,
1L, 2L, 1L, 2L), .Label = c("Ly49-CD8+TCell", "Ly49+CD8+TCell",
"SPtetramer+CD8+TCell"), class = "factor"), Ly = c("SPtetramerplus",
"SPtetramerplus", "SPtetramerplus", "SPtetramerplus", "SPtetramerplus",
"Ly49minus", "Ly49plus", "Ly49minus", "Ly49plus", "Ly49minus",
"Ly49plus"), treatment = structure(c(2L, 2L, 1L, 1L, 1L, 2L,
2L, 2L, 2L, 2L, 2L), .Label = c("MOG", "MOGSP"), class = "factor")), class = "data.frame", row.names = c(NA,
-11L))


Timeout when downloading the ncbi nr blast database I am experiencing timeout problems when downloading the NCBI nr preformatted blast database using the update_blastdb script (version 504861).
I run the script with the following paramters
update_blastdb --decompress --passive --verbose nr

and I get the following error message (in verbose mode)
Downloading nr (45 volumes) ...
Downloading nr.00.tar.gz...Net::FTP=GLOB(0x5610fb59b8f8)>>> PASV
Net::FTP=GLOB(0x5610fb59b8f8)<<< 227 Entering Passive Mode (165,112,9,229,195,144).
Net::FTP=GLOB(0x5610fb59b8f8)>>> RETR nr.00.tar.gz
Net::FTP=GLOB(0x5610fb59b8f8)<<< 150 Opening BINARY mode data connection for nr.00.tar.gz (18745730730 bytes)
Net::FTP: Net::Cmd::getline(): timeout at /usr/share/perl/5.26/Net/FTP/dataconn.pm line 82.
Unable to close datastream at /usr/bin/update_blastdb line 202.
Net::FTP=GLOB(0x5610fb59b8f8)>>> PASV
Net::FTP: Net::Cmd::getline(): unexpected EOF on command channel: Connection reset by peer at /usr/bin/update_blastdb line 203.
Failed to download nr.00.tar.gz.md5!
Net::FTP: Net::Cmd::_is_closed(): unexpected EOF on command channel: Connection reset by peer at /usr/bin/update_blastdb line 101.
Net::FTP: Net::Cmd::_is_closed(): unexpected EOF on command channel: Connection reset by peer at /usr/bin/update_blastdb line 101.

The timeout happens after ~35 minutes and a file that is approximately 18GB big is being downloaded, which matches the expected filesize. The checksum file (nr.00.tar.gz.md5) is not downloaded. So I'm not sure which of the two files is actually the problem.
I tested downloading the nt database and everything seems to work fine, so I don't think the script is the problem. For comparison, this is the nt download output for the first file
Downloading nt.00.tar.gz...Net::FTP=GLOB(0x562575a73168)>>> PASV
Net::FTP=GLOB(0x562575a73168)<<< 227 Entering Passive Mode (165,112,9,229,196,51).
Net::FTP=GLOB(0x562575a73168)>>> RETR nt.00.tar.gz
Net::FTP=GLOB(0x562575a73168)<<< 150 Opening BINARY mode data connection for nt.00.tar.gz (4065912989 bytes)
Net::FTP=GLOB(0x562575a73168)<<< 226 Transfer complete
Net::FTP=GLOB(0x562575a73168)>>> PASV
Net::FTP=GLOB(0x562575a73168)<<< 227 Entering Passive Mode (165,112,9,229,195,107).
Net::FTP=GLOB(0x562575a73168)>>> RETR nt.00.tar.gz.md5
Net::FTP=GLOB(0x562575a73168)<<< 150 Opening BINARY mode data connection for nt.00.tar.gz.md5 (47 bytes)
Net::FTP=GLOB(0x562575a73168)<<< 226 Transfer complete

Any help would be appreciated.

Same here. I used axel to download it a week ago, but fail to get finished. Also, the downloaded taxa zipped file (60Mb) doesn't match md5sum. Apparently, they update the FTP files every 30 minutes or 1 hour.


Gene Ranking - signal to noise ratio used in GSEA-P algorithm? I'm looking at Broad Institute's orignal GSEA-P algorithm R script which I downloaded here: http://software.broadinstitute.org/gsea/downloads.jsp.

I'm trying to adapt their GSEA.1.0.R script to process datasets that have 1 gene expression profile as for class 1, and N gene expression profiles for class 2. The bottleneck seems to be the GSEA.MetaboliteRanking() method. I looked and it appears they are using "signal to noise ratio" (SNR) metric to rank genes in patient profiles.

I was wondering how the signal to noise is calculated from gene expression profiles?

Note: I am using/modifying their script because it's the only implementation of the GSEA algorithm that I could find where they expose their code.

i.e., what is the mathematical reasoning behind the following lines of code? Note, there are little to no comments in the code, but the following variables, from what I can tell, are:

A: Matrix of gene expression values (rows are genes, columns are samples)
class.labels: Phenotype of class distinction of interest. A vector of
              binary labels having first the 1s and then the 0s

P = class.labels1 * subset.mask
n1 = sum(P[,1])
M1 = A %*% P
M1 = M1/n1
A2 = A*A
S1 = A2 %*% P
S1 = S1/n1 - M1*M1
S1 = sqrt(abs((n1/(n1-1)) * S1))

P = class.labels2 * subset.mask
n2 = sum(P[,1])
M2 = A %*% P
M2 = M2/n2
A2 = A*A
S2 = A2 %*% P
S2 = S2/n2 - M2*M2
S2 = sqrt(abs((n2/(n2-1)) * S2))

The noise is comming from the subset.mask, which is created above in a loop with the number of permutations.

for (r in 1:nperm) { #L90
....
        subset.mask[, r] <- as.numeric(c(subset.class1, subset.class2)) # L107

So by multiplying the random selection of subsets by the expression we get the "noise". Later on line 251 we get the ratio of noise/signal by dividing the signal from the contrast and the signal by the noise.

Those transformation on the middle seem to come from normalizing the signal in order to be able to compare the resampling, but I don't fully understand why they do this. Also note that the R scripts are outdated and no longer maintained, so they can contain errors.


TCGA bam file duplicates (aligned reads) and related last part of TCGA barcode I am running an analysis on TCGA bam files (GCRh37/hg19) of whole exome sequencing data.

I usually find one bam file for each sample. However, at least one sample has two bam files, in which the last part of the name (before the extension) has the suffix .1 or .2 or .3...

For example, C500.TCGA-BL-A3JM-01A-12D-A21A-08.1.bam and C500.TCGA-BL-A3JM-01A-12D-A21A-08.2.bam.

I didn't find an explanation for this last suffix. Are they duplicates, or.. ? Is there some wiki about the last part of the filename?

If you check the bam list of TCGA, you'll see that they use different numbers for the bam files (maybe different settings). What the difference is, I don't know, but if you want the file that they used for their paper, it is C500.TCGA-BL-A3JM-01A-12D-A21A-08.2.bam, the other one is not in the list.


STAR Indexing Diference for Small and Large Genome File Output I have a quick doubt on the output of the Genome Indexing, I have used the STAR program along with genome .fasta file and GFF file.

Genome size is 3GB, here is the file output

chrLength.txt
chrNameLength.txt
chrName.txt
chrStart.txt
genomeParameters.txt

I have another small Genome 60MB in size, I did the genome indexing, here is the file output

chrLength.txt
chrNameLength.txt
chrName.txt
chrStart.txt
exonGeTrInfo.tab
exonInfo.tab
geneInfo.tab
Genome
genomeParameters.txt
SA
SAindex
sjdbInfo.txt
sjdbList.fromGTF.out.tab
sjdbList.out.tab
transcriptInfo.tab

My point here is that, why I got the extra information for my small genome size, but I didn't get the same for the big size genome. I do apply the same procedure for the both.

here is the below information. Only difference I made for the large Genome size is (--sjdbOverhang 99 --genomeChrBinNbits 15) to reduce the memory, but the rest of things are same for small genome.

#!/bin/bash
NUM_THREADS=12
mkdir DB
STAR --runMode genomeGenerate --genomeDir DB \
    --runThreadN $NUM_THREADS \
    --genomeFastaFiles XL9_2.fa \
    --sjdbGTFfile XENLA_Frog.gtf+gff3 \
    --sjdbOverhang 99 \
    --genomeChrBinNbits 15

Could anyone give an idea, why there is different, I am new to this field, so I am wondering about the difference in this.

Thanks in advance.

Cheer
San

You didn't get extra information with the small genome, the index for the larger one just wasn't made. My guess is that you ran out of RAM. For mammalian sized genomes ensure you have at least 30GB of RAM, likely more like 60GB. Also, make sure to look for error messages.


Lipidomics normalization (peak area) I have data from a lipidomic experiment comprising: 2 conditions and 5 samples in each condition; and for each lipid the value is a peak area.
I am doing a lipidomics analysis using the web app lipid ontology, and they recommend to normalize the data prior to doing the analysis,

"There are many ways to normalize data, so LION/web cannot check
whether data is properly normalized. If you don't know how to
normalize the data, we recommend normalization by expressing all lipid
species as fraction of the sum. A popular website that supports
several normalization algorithms is MetaboAnalyst."

I don't know how to normalize my data, I cant find any paper related to lipidomics normalization to know to best way to do that, and also I can't find normalization function in the MetaboAnalyst.
Thank you in advance!

Basically and briefly,

Normalisation = transforming a data set to conform to a normal distribution. Converting to a frequency is NEVER a normalised transformation
Standardisation = ensuring each point of the data set is comparable to the others, either directly or relatively depending on the calculation

The information you have received is incorrect, moreover you only convert to a frequency for a final assessment of the data and in my experience NEVER as an intermediate step. Normally data point between 0 and 1 are bad news and often require a transformation to get them away from fractions.
I don't understand the behaviour of the variance involved in metabolism and particularly lipids, however you should refer this issue to "cross-validation" Stackexchange. I did have some understanding of enzyme kinetics but in a different area.
Cross-Validation Stackexchance statisticians will understand less of the biology than I do (which is a bit worrying of course), so you will need to give them quite detailed information on the data set. They will almost certain support my concern of the advice you are given and provide a better method of standardisation. One method is to take the difference in mean and divide it by the standard deviation, you can then normalise the standardised data. They will want to know the goals of your investigation to understand who to process the data further.
What I can say is variation between runs is often an issue (in your apparatus perhaps not, I dunno) and standardising against a common control is one approach to overcome this, e.g. the control same common to all runs would be the mean about which all other data is centred.
What I suspect is that lipids are less popular than other fields, so standard protocol may not be observed.


Difference between 1D and 1D^2 data I am working on ONT data. Initially, I have worked on data from 1D using Minimap2 aligner. I came across 1D^2 from ONT website. I was wondering how different is the data from two techniques? What is the difference between 1D and 2D? and do we get 2  fastq files (one per strand) for each sample?

TIA

I would greatly appreciate any help regarding my question.

1D² data that is called by Albacore generates two overlapping sets of files. Here's what Devon O'Rourke mentioned on the ONT community forums and on Twitter:

  Following Albacore basecalls on a 1d2 library I get two sets of .fq files, summary stats, etc.: one for the 1D basecalling script, and one for the 1D2 script. I'd like to use the reads generated from these scripts in a genome assembly.

  After a little bit of grep searching, it seems like there are overlapping reads among the 1D and 1D2 directories - namely, those 1d reads that generated 1d2 reads!

Link here, for those who have ONT community access.

I've recommended generating a list of the 1D^2 reads, then excluding those reads from the 1D fastq files. I have a couple of scripts that can help for this:

$ cat 1d2_reads.fq | ~/scripts/fastx-length.pl | perl -lane '$F[1] =~ s/(.{36})/$1\n/; print $F[1];' > readNames.txt
$ cat 1d_reads.fq | ~/scripts/fastx-fetch.pl -v -i readNames.txt > filtered_1d_reads.fq

More information on this can also be found in this StackExchange answer.


Reading fasta sequence regions with faidx Im coding in c++ and reading in different reference genomes to examine regions across the chromosomes a few hundred basepairs at a time.
To do it in c++ i use the library htslib and the command faidx_fetch_seq() which is similar to read specific regions using samtools faidx.
The problem is when im reading in the reference genomes, there are of course multiple regions consisting only of NNNNN (and so on). Since i only need the actual DNA sequence i wish to read in the entire file but removing the regions with NNNN.
I have tried reading in the regions in c++ and then just creating an if-statement to check the actual bases consist of N. This works but is incredibly slow.
So my question is: do anyone know how to only read in the actual DNA sequences, removing the NNN regions? either with samtools or with the c++ funciton faidx_fetch_seq().

If you can't read in the entire sequence for a given reference genome:
(1) Use a script to generate a BED file of masked regions of low complexity (i.e. where you find Ns in your FASTA file). (2) Use bedops --complement with this mask file and the output of UCSC fetchChromSizes to generate a BED file containing regions that are unmasked. (3) Repeat for other reference genomes, as needed.
You can now loop over the unmasked regions in this file, making calls to samtools faidx or faidx_fetch_seq(). These calls will return unmasked sequence.
For performance or convenience, you could read chunks of sequence into a std::deque or std::vector buffer and then read through that however many hundreds of bases, refreshing the buffer, as needed.


What is the standard way to measure contig sequence lengths in a BAM? What is the standard way to measure contig sequence lengths in a BAM?

My understanding is that the community would use samtools idxstats to compute this information from the corresponding index file.

Are there more precise/standard approaches?

Standard approaches would be either samtools idxstats or samtools view -H, both of which will produce the exact same results. In fact, all methods will produce the exact same result, since contig length is set by the reference you align against.


.SDF file to .PDB file I am currently working on identifying effective antiviral compounds against the corona virus. I have downloaded the Covid-19 AntiViral Candidate Compounds Dataset from the www.cas.org webpage, the link to which I have provided below. Me and a few other doctors are using Docking software to predict the docking capacity of these anti-viral compounds. Docking software requires the anti-viral compounds to be in the .pdb file format only then is it recognized by Automated docking software. Is there any way I can convert the .sdf file containing the anti-viral compounds to .pdb file or is there a way of obtaining the .pdb files for list of anti-viral compounds.

Thank you for your time and patience,
Anas Hussain.

Anti-Viral Candidate Compounds Dataset
https://www.cas.org/covid-19-antiviral-compounds-dataset

Use OpenBabel.
Something like babel -isdf file_in.sdf -opdb file_out.pdb.


Clustering information saved in Seurat object Seurat V2 had a option to find clustering information saved in object: PrintFindClustersParams(object = pbmc). How can I get the same clustering parameters from objects in Seurat3? Do I need to manually save clustering parameters or they are saved automatically?
Update: the command line object@commands works perfectly for ...$FindClusters. However, if I would like to find information about dimensions used there is an error:
...$FindNeighbors.SCT.pca Error in [[.Seurat(x, i, drop = TRUE) : Cannot find 'FindNeighbors.SCT.pca' in this Seurat object

I am not sure of the output of PrintFindClustersParams, as I haven't use it before, but try with object@commands. It returns a list containing every command and parameters used.
For example
...
$FindNeighbors.SCT.pca
Command: FindNeighbors(seurObj, dims = 1:20, verbose = vrb)
Time: 2019-12-11 12:41:12
reduction : pca
dims : 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
assay : SCT
k.param : 20
compute.SNN : TRUE
prune.SNN : 0.06666667
nn.method : rann
annoy.metric : euclidean
nn.eps : 0
verbose : FALSE
force.recalc : FALSE
do.plot : FALSE
graph.name : SCT_nn SCT_snn

$FindClusters
Command: FindClusters(seurObj, verbose = vrb)
Time: 2019-12-11 12:41:17
graph.name : SCT_snn
modularity.fxn : 1
resolution : 0.8
algorithm : 1
n.start : 10
n.iter : 10
random.seed : 0
group.singletons : TRUE
verbose : FALSE
...

Re: update
I think it depends on how you built the object. FindNeighbors.SCT.pca should be result of running sctransform normalization. If you use the standard way, you may have other names (e.g.FindNeighbors.pca).


RNA-seq analysis of mixed viral/host reads with salmon I did RNA-seq of mammalian cells infected with pox virus. Now, I have read files which contain both host and virus reads. I want to align the reads both to host and viral genome. I was thinking I could concatenate the host and virus genome into one file and run salmon against this concatenated genome. However, salmon recommends transcriptome file for assembly which are not available for viruses. Virus genome are available as GenBank or GFF3 format from NCBI. Is there any way I can concatenate these formats into the format that can be used by salmon? Or is there any way around to use virus genome as reference in salmon?

You provide Salmon with a transcriptome fasta... so merging the human transcriptome with the pox virus genome fasta file should work. You don't mention that in your post but the viral genome should exist as a fasta. You can extract the sequences from genbank or maybe you can find our virus here: https://www.ebi.ac.uk/genomes/virus.html

You can find the coding transcriptome at the link below in the cdna folder.
The non-coding transcriptome is in the ncrna. You should probably merge them:

ftp://ftp.ensembl.org/pub/release-95/fasta/homo_sapiens/


How to compare the contents of a column of the same data table I have this table:

and I want to get the rows that are equal to the first three columns, like this:

I've tried these functions, but when I get the index of the lines, r doesn't give the output that I want:

df$obj<-sapply(c("sample1", "sample2", "sample3"), function (comparison) {
    df$sample1 == df[comparison]
})

Loop through 3 columns by row using apply, then convert toupper case (as we want a same as A), get unique, and get lengths, if length is 1 then all 3 values are equal:

df1[ lengths(
  apply(df1[, 1:3], 1, function(i) unique(toupper(i)))
  ) == 1, ]

This would scale better if we want to compare 4, 10, n(?) columns.


DMPfold: Understand the run of DMPfold "run_dmpfold.sh"? I ran dmpfold using this command:

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

It produced files like this:

sequence_sub_embed_1.pdb ..... sequence_sub_embed_500.pdb

And, now I see some more files like this being produced:

sequence_1.pdb .... sequence_2.pdb .....etc

Can you please let me know what is dmpfold doing here? And can I stop dmpfold in between when desired number of models are generated? And, while it is running, can I move the previously made pdbs to new location becuase of disk space issues or all pdbs have to be in the "dmpfold-lysin" as stated in the command-line above? Thanks.

From the docs:

  The final model is final_1.pdb and other structures may or may not be generated as final_2.pdb to final_5.pdb if they are significantly different. Many other files are generated totalling around 100 MB - these should be deleted to save disk space if you are running DMPfold on many sequences.

These other files include intermediate CNS files (sequence_sub_embed_1.pdb etc.), output files before clustering (sequence_1.pdb etc.) and saved ensembles of previous iterations (ensemble.1.pdb etc.).

Let the run finish or you won't get the final output files. If you want to use all the generated structures, use the ensemble.n.pdb file from the last iteration (i.e. the highest number). You can also use the other ensemble files if you want to include models from previous iterations.

I would strongly recommend using default "3 50" parameters until you get things working. Don't move files during the run. You really don't get much benefit of going beyond 50 models per iteration, trust me. We tried it relatively extensively!


Is it important to filter out poor quality cells before performing an integration analysis on single cell RNA sequencing data? In order to perform an integration analysis of single cell RNA seq data, is it important to check the percentage expression of mitochondrial genes of cells as well as the feature counts to exclude dying or dead cells as well as multiplets from further analyses?

I ask because the standard Seurat work flow to perform integration as outlined in Integration and Label Transfer does not mention filtering steps prior to integration.

In this publication, Eleanor Alenzi et al. 2018 also did a canonical correlation analysis with Seurat but did not mention steps to filter out poor quality cells before doing the integrated analysis. This second publication does same. Both publications do not mention if their visualisations or differential gene expression analyses are done on integrated (discouraged) or raw uncorrected RNA data (recommended).

Following the styles of the publications linked above (many more do same), is it also efficient to perform integration without filtering out poor quality cells?

It is absolutely necessary to remove low quality cells: In the case of CCA (and this applies to other "integration" or "data alignment" methods as well), one would need to use "anchors", basically same type/state/kind of cells from the samples to be "integrated" and are used to "align" the different samples. In a scenario where you have a set of cells whose distinguished feature is high mitochondrial content due to apoptosis irrespective of their origin (epithelial, T cell, ...) and such cells are picked up as so called anchors, the whole integration effort will be off. Even when not using integration methods, bad quality cells quite often form their own clusters (highly skewed biology) not contributing anything to analyses and are often removed at the downstream steps anyway.

And when it comes to published work, unfortunately it comes down to the backgrounds of the reviewers. I tend to think that there are way too many single cell studies than experts with required bioinformatics and statistics so that a thorough review can be conducted. And lack of a paragraph on QC does not mean no QC has been performed anyway.


How to model the effect of acid/base on pH I would like to model a biological compartment with sysBio R package. It is fine but I'm having a difficult time with modeling the change in pH due to acidosis. I tried several things like:

addMAreaction(glu, "pH = null", "acidosis", "base excess")
addMAreactRate(glu, "acidosis", "assigned", "ka*pH")
addMAreactRate(glu, "base excess", "assigned", "ka*pH*-1")

which caused continuously decreasing pH levels. I'll leave out the intermediate steps but the simplified version of the last code looks like this:

library("sysBio")

glu <- newModel("glu")
addSpecies(glu, "Lactate", 0)
addSpecies(glu, "Glucose", 1000)
addSpecies(glu, "pH", 7.4)
addSpecies(glu, "Acid", 0)

addMAreaction(glu, "Glucose = Lactate + 0.01*Acid", "forwardL",
              "backwardL", name = "glycolysis, sort of")

addMAreactRate(glu, "forwardL", "assigned", "kgl*Glucose")
addMAreactRate(glu, "backwardL", "assigned", "ka*pH+klg*Lactate")

addParameters(glu, "kgl", 0.01)
addParameters(glu, "klg", 0.01)
addParameters(glu, "ka", 0.05)
addParameters(glu, "a", 0.001)

addRule(glu, "pH rule", "ODEs", "pH=7.4-a*Acid")

makeModel(glu)
simResults <- simulateModel(glu, times = seq(0,10,0.1))
plotResults(simResults, title="glu")

And the resulting plot is this where pH is increasing forever despite Acid is at a steady state concentration.

Please note that the abovementioned code is just a simplification. It doesn't even look like aerobic glycolysis where 1 glucose produces  lactate and 2 hydrogen atoms etc.

I intend to include several acid- and base-creating reactions and evaluate the resulting pH. So could someone point me in the right direction about how to model the relationship between pH and acids or bases?

The answer does not need to be explicit code. I'll appreciate any solid comment on my design errors.

I figured that one has to include only the changes in the rules formula. Keeping the initial value of the Species or Parameter when defining circular or self-referencing formulas caused the error. So the correct rule becomes:

pH=-a*Acid instead of pH=pH-a*Acid or pH=7.4-a*Acid.

For convenience, I include the correct code crowded with several reactions producing both acids and bases, where pH gradually decreases as expected. The names are fictional, and just are better placeholders instead of x1 and x2, etc.

library("sysBio")

glu <- newModel("glu")
addSpecies(glu, "Glucose", 1000)
addSpecies(glu, "Lactate", 0)
addSpecies(glu, "Pyruvate", 0)
addSpecies(glu, "Oxalate", 0)
addSpecies(glu, "Acid", 0)
addSpecies(glu, "Base", 0)
addSpecies(glu, "pH", 7.4)

addMAreaction(glu, "Glucose = Lactate + Acid", "forwardL",
              "backwardL", name = "glycolysis")
addMAreaction(glu, "Lactate = Pyruvate + Acid", "forwardP",
              "backwardP", name = "lactolysis")
addMAreaction(glu, "Pyruvate = Oxalate + Base", "forwardO",
              "backwardO", name = "oxalate synthesis")

addMAreactRate(glu, "forwardL", "assigned",
               "k_glucose_to_lactate*Glucose")
addMAreactRate(glu, "backwardL", "assigned",
               "k_acid*pH+k_lactate_to_glucose*Lactate")
addMAreactRate(glu, "forwardP", "assigned",
               "k_lactate_to_pyruvate*Lactate")
addMAreactRate(glu, "backwardP", "assigned",
               "k_acid*pH+k_pyruvate_to_lactate*Pyruvate")
addMAreactRate(glu, "forwardO", "assigned",
               "k_pyruvate_to_oxalate*Pyruvate")
addMAreactRate(glu, "backwardO", "assigned",
               "k_base*pH+k_oxalate_to_pyruvate*Oxalate")

addParameters(glu, "k_glucose_to_lactate", 0.02)
addParameters(glu, "k_lactate_to_glucose", 0.001)
addParameters(glu, "k_lactate_to_pyruvate", 0.01)
addParameters(glu, "k_pyruvate_to_lactate", 0.0001)
addParameters(glu, "k_pyruvate_to_oxalate", 0.01)
addParameters(glu, "k_oxalate_to_pyruvate", 0.001)
addParameters(glu, "k_acid", 0.0001)
addParameters(glu, "k_base", 0.0001)

addRule(glu, "pH rule", "ODEs", "pH=-k_acid*Acid+k_base*Base")

makeModel(glu)
simResults <- simulateModel(glu, seq(0,3,0.01))
plotResults(simResults, title="glu")

The head of the resulting data.frame is:

  time   Glucose  Lactate   Pyruvate      Oxalate     Acid         Base       pH
1 0.00 1000.0000   0.0000  0.0000000 0.000000e+00   0.0000 0.000000e+00 7.400000
2 0.01  844.0001 155.0615  0.9384259 1.402282e-05 156.9384 1.402282e-05 7.399914
3 0.02  785.0869 210.3646  4.5476263 7.892761e-04 219.4615 7.892761e-04 7.399719
4 0.03  772.1557 218.6462  9.1924839 5.655403e-03 237.0425 5.655403e-03 7.399489
5 0.04  768.8212 217.2720 13.8876293 1.917807e-02 245.0856 1.917807e-02 7.399248
6 0.05  767.0383 214.4883 18.4278053 4.550630e-02 251.4350 4.550630e-02 7.399000

and resulting graph looks like:


Building network for biomarker interaction analysis I am looking for advice on how to best approach a problem I am faced with.

I have a dataset of numerous degradative biomarkers, clinical information and various other measures (from clinical trials). I would like to see how different biomarkers relate to each other, and effect their levels (high/low etc).

My goal is to gain an understanding of how these biomarker levels fluctuate in different disease profiles - and understand if there is a relationship between each of them (co/independent). The limitation is that I do not have access to genetic or environmental data so this will be a "random effect". There are around 1000 marker samples and 10's of biomarkers mixed between clinical and biochemical. There are also numerous time points to be incorporated.

I have thought about implementing something similar to a protein-protein interaction network, although this will of course not involve physical interactions, but more pathway interactions for further hypothesis driven research.

I wanted to firstly: validate if this project actually makes sense and secondly: ask for some advice around how to best implement this.

I have experience implementing neural/deep nets, directed and undirected networks (Bayes) in python, but not sure exactly how to go about this.

Thanks!

You can start by looking at the correlation (if data is non-parametric) of each variable with each variable of each data point you have. One thing that might change is the correlation between variables along the disease.

Next you can determine which relationships are constant along the disease progression/clinical trial. Identify also the ones that don't show a constant relationship. To do so you can use WGCNA (although with just 10 variables is too low to expect a scale free network you can use to correlate the variables and group them by correlation profile)

Then the question is which variables are more important, those who kept the same relationship or those which change?
Use the eigengenes of those variables (separately by those that have a constant relationship and those which don't) to estimate the dependency with the parametric variables. Or you simply can do an ANOVA of all the variables with each non-parametric variable.


What does liability mean in GWAS heritability? I am reading about GWAS in heritability. They usually say to calculate heritability on a liability scale but after searching for this word "liability", I still don't understand clearly what does liability mean? Does anyone have any easy-to-understand explanation regarding this term? Thank you very much!

When you use linear mixed models to estimate heritability you assume that the underlying trait is normally distributed which is called a disease liability scale.

For continuous traits this is not a problem but for binary traits, this becomes an issue because you have a 0/1 value for a phenotype and usually there is a higher proportion of cases in the study sample than the general population prevalence of disease, which leads to an ascertainment bias.

The population prevalence also varies between populations, for example the prevalence of malaria in one continent is vastly different from another.

So to make the heritability estimates comparable and also to correct for the ascertainment bias, the observed scale heritabilities for dichotomous traits are transformed to the liability scale taking population prevalence and the case proportion into account.

Some good papers for reading:

Falconer D.S. The inheritance of liability to certain diseases, estimated from the incidence among relatives. Ann. Hum. Genet. 1965;29:51–71
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3059431/


How to map short sequences to long reads, recovering all multiply-mapped high-quality matches The dilemma:
I have around one million short sequences (21 bp to several 100s of basepairs) for which I need to identify all occurrences of in 20-30x coverage noisy long reads (both pacbio and ONT). All of the short sequences and long reads are derived from the same individual, and there are enough short sequences given the organism's genome size that each long read should have at least one or multiple short sequences.
An example:
Given three short sequences, AA & BBBB & CCCCCC , how to identify all occurrences of these sequences in long, uncorrected reads.
Sequences to identify in long reads:
  - AA
  - BBBB
  - CCCCCC

Long Reads:
  - 1: ------AA-----------------------CCCCCC------
  - 2: ------AA--------BBBB-----------------------
  - 3: ---BBBB-----------CCCC--
  - 4: ------------------AA-----

In this example, the output sam file will contain high MAPQ hits for short sequence AA to reads (1, 2, 3), for short sequence BBBB to reads (2, 3). For short sequence CCCCCC there will be a high MAPQ hit to read 1 and a lower MAPQ hit to read 3
Limitations:

Due to the requirements of this project, I cannot perform a de novo assembly then map the short sequences to the assembly output. It is necessary to directly map the short sequences to the long reads.
I cannot correct the long reads before I map/align the short sequences to them.
The alignment/mapping should sensitively detect the short sequences in the long reads.
The short sequences may be anywhere from 1Mbp to 30Mbp in total bases.
The long reads may be anywhere from 3Gbp to 90Gbp in total bases. So the mapper/alignment technique is hopefully fast. Sensitivity is more important though.

Research so far:

bwa mem does not output multiply-mapped in a predictable manner that fits this use case https://www.biostars.org/p/304614/
blat seems to output many hits for a single short sequence, but there are no out-of-the-box parameters for mapping short reads to long reads.
Maybe multiple sequence alignment is the most sensitive, but wouldn't that entail running (no. of short seqs) * (no. of long reads) alignments?

LAST has given the best results for me when I've tried to do this, although I agree with @user172818 that it's not a good idea to map really short reads. This is due to a combination of natural sequence duplication in long DNA sequences (e.g. see here), as well as abundant base calling differences present in single-molecule sequencing.

Minimising error is not necessarily going to be the best option, and concentrating only on unique hits will miss a lot of signal. There are frequently multiple identically plausible positions, even at zero error, and the long reads will have their own associated errors.

I also find the limitations a little odd. Canu can do read-level correction on long reads based on overlapping reads, and if assembly is possible, then nanopolish can be used to correct nanopore reads for systematic base-calling error introduced by the methylation of unamplified template.


How to query NCBI (Nucleotide database) by a feature qualifier? I cannot seem to find a way to query nucleotide database (https://www.ncbi.nlm.nih.gov/nuccore) by a specific qualifier of the feature. For example, by /host under the feature "source".
For example, below is a part of one record, where the qualifier /host contains "Bemisia tabaci B".
How could I search for all sequences that contain the word "Bemisia" in the qualifier \host?

FEATURES             Location/Qualifiers
 source          1..351658
                 /organism="Candidatus Portiera aleyrodidarum BT-B-HRs"
                 /mol_type="genomic DNA"
                 /strain="BT-B"
                 /host="Bemisia tabaci B"
                 /db_xref="taxon:1206109"
 gene            1..1926
                 /gene="dnaK"
                 /locus_tag="C530_RS00005"
                 /old_locus_tag="C530_001"
 CDS             1..1926
                 /gene="dnaK"
                 /locus_tag="C530_RS00005"
                 /old_locus_tag="C530_001"
                 /inference="COORDINATES: similar to AA
                 sequence:RefSeq:WP_014894920.1"
                 /note="Derived by automated computational analysis using
                 gene prediction method: Protein Homology."
                 /codon_start=1
                 /transl_table=11
                 /product="molecular chaperone DnaK"
                 /protein_id="WP_014894920.1"
                 /translation="MIRIIGIDLGTTNSCVAVLENGKSRVIENSEGARTTPSIIAYTD
                 DNEILIGQTAKRQAVTNPSNTLYAIKRLIGRKFNDDIVQKDIKMVPYKIVSAENGDAW
                 VQVKNRKLAPPQISAEILKKMKKTAEDYLGEKITEAVITVPAYFNDSQRQATKDAGRI
                 AGLNVKRIINEPTAAALAYGMDKNRGDKKIAVYDLGGGTFDISIIEIADVDGETQFEV
                 MATNGDTFLGGEDFDLKLINYIVKQFKLESGVDLSGDSLAMQRLKEAAEKGKIELSSS
                 QQTDINLPYITADKTGPKHLTLKITRAKLESLVEDLIKRSLLPCKKAIDDAKLTTKEI
                 DDIILVGGQTRMPLVQKKVAEFFGKEARKDVNPDEAVAMGAAIQGGVLGGDVKDVLLL
                 DVTPLTLGIETMGGIMTPLIEKNTTIPTKKTQIFSTAEDNQTAVTIHALQGERKKALQ
                 NKSLGRFDLTDIPPAPRGTPQIEVAFDLDANGILNITAKDKATGKEQSIIIKSSGGLS
                 EEEINKMIQDAETNKEEDKKFEELVQTRNQADGMIHTIRKTLKEIKNEDEKNKLESLI
                 NNLEQTLKSDDIKTIKEDLNKLTEEYSTIYQKIYSEKKPQNHNSEPAQKNAEKSSSKK
                 EGDVVDAEYEEVNNKTK"
                 ...

Best,
Milan

For this particular query, try host=Bemisia tabaci B[WORD].  However, you may also be interested in related queries like host=Bemisia tabaci B biotype[WORD].
But, I am not sure a general solution exists, as I do not believe all qualifiers are indexed.

Also, you can see how your query is interpreted on the right side under "Search details".  For the original attempt ([host]"Bemisia tabaci B"), it is interpreted as:
host[All Fields] AND "Bemisia tabaci B"[All Fields]

while the goal is to produce
host=Bemisia tabaci B biotype[WORD]


Processed spliced and unspliced count matrices for existing scRNA-seq atlases Before I do this myself...

Is there any compiled batch-corrected dataset of the major scRNA-seq atlases (Mouse Cell Atlas, Mouse Organogenesis Atlas, Mouse Gastrulation Atlas, Tabula Muris, Tabula Muris Senis, etc.) with processed spliced AND unspliced count matrices?

I would be surprised if there is, but in the event there isn't, are there new pipelines in development that will make processing of raw data easier for studying RNA velocity?

I am not sure if there are preprocessed datasets available, however, Kallisto + Bustools was designed specifically for the rapid quantification of spliced and unspliced transcripts.


How can I match sequences with IDs using python dictionaries? 0.980634         -                    572         sp|P06672|PDC_ZYMMO                         -                     568                               1.2e-263             874.5          0.0            1        1      6.6e-267  1.4e-263   874.2      0.0       5         562    15     570    12     572    0.98   Uncharacterized                         protein                   OS=Lomentospora            prolificans        OX=41688             GN=jhhlp_005678  PE=3             SV=1

I've asked a similar question before. I have a better understanding of python now and would like to revisit this. Can someone critique my thinking?
I have an HMM output file that contains about 20 fields. There are two sets of IDs, e.g., P006672 and jhhlp_005678
In a separate file are FASTA-formatted sequences linked to these IDs.
Proposed solution:
1
   a) read in fields containing IDs, lets called them fields 8 and 6;

   b) store field values in an array or dictionary (whatever the proper data structure is for
  this). Let's called this Fields_8_and_6

a) read proteome file containing FASTA-formatted sequences
b) Grep and store as ID:sequence key:value pairs; let's call this ID_to_sequence

a) Compare Fields_8_and_6 to ID_to_sequence
b) Write ID1:sequence1 ID2:sequence2 as pairs into one file

send to ClustalO

Biopython can index a fastafile and you can access the sequences by their ID in a dict-like manner.
All that is left for you to do is extract the IDs from your file. (more Biopython)


Pooling data in metagenome assembly I have 12 human gut microbiome WGS Nextseq reads (151 bp paired end). What will be an effective strategy to assemble a metagenome?

Let us say I have already filtered the fastq for quality, adapter sequence and host contamination (human, in this case).

1) Should I concatenate all the R1 reads as one single R1 read and one Single R2 read?

cat Sample[1..12].R1 > Single_R1.fastq
cat Sample[1..12].R2 > Single_R2.fastq

and then use Diginorm to normalize the Single_R1.fastq and Single_R2.fastq. Subsequently, feed these fastq files into any metagenome assembler such as Megahit, MetaSPAdes?

Normalize the output by using CD-HIT or similar tool to remove duplicates and filter by contig length.

OR

2) Perform metagenome assembly for each of the samples individually after applying filtering, removing adapters and host contamination.

R1=(*_R1_001.filtered.fastq)
R2=(*_R2_001.filtered.fastq)
for ((i=0;i<=${#R1[@]};i++)); do
  /bin/metagenome-assembler -1 "${R1[i]}" -2 "${R2[i]}" -o ${R1[i]%.*}.contigs.fa;
done

Followed by combining all the contigs.fa into one mega_contigs.fa

cat *.contigs.fa > Mega_contigs.fa

and use CD-HIT or similar tool to remove duplicates.

With metagenomes, there's several strategies and no one-size-fits-all.
 You can spend a lot of time playing with the data to try to get the "best assembly."  But it depends on what you're trying to do, e.g. whether or not you are mostly interested in the highest abundant species or are interested in the less abundant ones too.

But to propose a strategy that seems reasonable and hasn't been mentioned yet, my initial thoughts are that since the gut microbiome isn't complex, you can likely get decent assemblies from each run separately.  Then you can run a contig clustering algorithm (e.g. MetaBAT) to bin the contigs into genomes.  Next compare the bins from the various assemblies (i.e. cluster the clusters).

It may or may not be worth the trouble to try to combine the related bins' contigs into a consensus (pan)genome sequence.  If you want to, methods that have been used successfully include shredding the contigs into overlapping chunks and assembling with an OLC assembler.  But if you're interested in KO:read-counts (for example), a consensus genome sequence isn't really required.  Cross-mapping the unassembled reads can be useful too (i.e. a sample may have a low-abundance OTU that doesn't produce good contigs, but it can be associated with an OTU that is better assembled in another sample).

If you provide more detail on your goal, I can comment further.


Read alignment using Bowtie2 So this is related to CRISPR-CAS9. I am working with off-target predictions for my thesis and was looking at all scientific papers related to CRISPR. I found one and decided to use their datasets. The datasets were not available directly and the method used to obtain them was specified which I had mentioned below:

They have around 30 sgRNAs which are each 23bp long and are aligning the sgRNAs to human reference genome with a tolerance of upto 6 mismatches using Bowtie2.

In order to access the exact dataset, I thought I need to replicate the method they used. To do so, I collected the sgRNAs from their paper and tried to use Bowtie2 to align them against human genome with upto 6 mismatches.

Since my laptop does not have the capacity to perform this operation, I am searching for interfaces that provide access to Bowtie2. I found Galaxy to be useful. Does anyone know of any other tool that provides an option for Bowtie2?

Thanks!

One possible workaround would be to blast your sequences against the genome. Then, when you have found where they align, you can extract those regions only into a new fasta file (make sure to take a few kb around each target region) and then use that fasta file as the "genome" you pass to bowtie2.

That would give bowtie a far smaller file to work with and should reduce the memory needed significantly.


Where to download baseline/average gene expression level of all human coding genes? I am looking for the most appropriate dataset for downloading baseline gene expression level across all human coding genes during development. I am aware that EMBL Expression Atlas is one of the resources that provide such information, but I am ideally looking for species averaged (and probably scaled) expression data for all human coding genes during early stages of development!
Any insight is much appreciated!

There is not going to be a public database with this data. Apart from anything else, generating data for early human development is difficult and ethically tricky. Also, most people who care about development would probably think that a whole body average would be meaningless.

The best you could do would be to select tissues you were interested in from something like Expression Atlas and calculate the averages yourself.


BUSCO does not find buscos in reference sequence of C. elegans (technical problem) I am working on a GNU/Linux machine. I have installed BUSCO using conda. I have downloaded the nematode dataset nematoda_odb9.tar.gzhere. I downloaded the reference sequence of C. elegans from Ensembl. I then ran busco using the following command:

run_busco --in Caenorhabditis_elegans.WBcel235.dna.toplevel.fa --cpu 1 --out busco --lineage_path ./nematoda_odb9/ --mode genome --force > busco.o 2> busco.e

tblastn seems to have run properly. run_busco/blast_output/tblastn_busco.tsv contains lots of hits for many of the buscos.

However, I have 982 "Missing BUSCOs".

The file that contains logs and errors related to augustus (run_busco/augustus_output/augustus.log) contains:

Will create parameters for a EUKARYOTIC species!
creating directory ...anaconda/envs/busco3_0_2/config/species/BUSCO_busco_1395218643/ ...
creating ...anaconda/envs/busco3_0_2/config/species/BUSCO_busco_1395218643/BUSCO_busco_1395218643_parameters.cfg ...
creating ...anaconda/envs/busco3_0_2/config/species/BUSCO_busco_1395218643/BUSCO_busco_1395218643_weightmatrix.txt ...
creating ...anaconda/envs/busco3_0_2/config/species/BUSCO_busco_1395218643/BUSCO_busco_1395218643_metapars.cfg ...
The necessary files for training BUSCO_busco_1395218643 have been created.
Now, either run etraining or optimize_parameters.pl with --species=BUSCO_busco_1395218643.
etraining quickly estimates the parameters from a file with training genes.
optimize_augustus.pl alternates running etraining and augustus to find optimal metaparameters.

It does not seem to contain any errors.

I have the following error:

Failed to open sequence file run_busco/augustus_output/extracted_proteins/EOG091H07ZC.faa.1 In fact, the folder run_busco/augustus_output/extracted_proteins/ is empty.

Anyone would know why I have no extracted proteins?

Most likely your Augustus installation is not functioning properly. Did you install Augustus from source or using conda? I believe the conda version of Augustus is known to have some issues.

Something worth trying: Running a complete busco install in a VM. This image has all dependencies installed and properly configured.


Chromatin accessibility level on promoters with different CpG ratio calculation I want to categorize promoters based upon high, intermediate and low-CpG content (high-CpG-density promoters (HCPs), intermediate-CpG-density promoters(ICPs), and low-CpG-density promoters (LCPs)).
So the data I have is for promoter which i have annotated 1000+/- around TSS and taken them as something less than 1kb is promoter region and beyond distal.
In term of tool i have used chipseeker to annotate and do the above step.
So now if i have to find the CpG density as i have mentioned above how do i get that?

Below I import the peaks call in narrowbed format from macs2 into a GRanges object. If you have a bed file, you can just use rtracklayer for it:
library(rtracklayer)
library(BSgenome.Hsapiens.UCSC.hg38)

peaks_narrowbed = "https://www.encodeproject.org/files/ENCFF846JAO/@@download/ENCFF846JAO.bed.gz"

extraCols <- c(signalValue = "numeric", pValue = "numeric",
                          qValue = "numeric", peak = "integer")
gr <- import(peaks_narrowbed, format = "BED",extraCols = extraCols)

head(gr)

GRanges object with 6 ranges and 6 metadata columns:
      seqnames        ranges strand |
         <Rle>     <IRanges>  <Rle> |
  [1]     chr1   10049-10326      * |
  [2]     chr1 180686-181040      * |
  [3]     chr1 186648-186903      * |
  [4]     chr1 191293-191604      * |
  [5]     chr1 267872-268146      * |
  [6]     chr1 586110-586324      * |
                                                   name     score signalValue
                                            <character> <numeric>   <numeric>
  [1] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_1       135     7.10907
  [2] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_2       309    11.65888
  [3] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_3        81     5.40289
  [4] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_4        56     4.54981
  [5] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_5       145     7.39343
  [6] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_6       116     6.54035

Then we extract the sequences that are within these regions, you can of course redefine the region by extending it from the peak:
library(Biostrings)
SEQ = getSeq(Hsapiens,gr)
head(SEQ)

A DNAStringSet instance of length 6
    width seq
[1]   278 TAACCCTAACCCTAACCCTAACCCTAACCCTAAC...CCCCAACCCCAACCCCAACCCCAACCCCAACCC
[2]   355 CTAACCCCTAACCCCTAACCCTAACCCTACCCTA...CTCAGCCGGCCCGCCCGCCCGGGTCTGACCTGA
[3]   256 GCAACCACCTGAGCGCGGGCATCCTGTGTGCAGA...CCCCTTCTTTCCATTGGTTTAATTAGGAACGGG
[4]   312 CAATCAGCAGGGACCGTGCACTCTCTTGGAGCCA...CCCTGCCCCTGTCTCTTCCGTGCAGGAGGAGCA
[5]   275 TGGATGTGTGCATTTTCCTGAGAGGAAAGCTTTC...TTGTTTAGTTTGCGTTGTGTTTCTCCAACTTTG
[6]   215 GAGGCTGAAGGAGACTGATGTGGTTTCTCCTCAG...CCAACTTTGTGCCTCATCAGGAAAAGCTTTGGA

This is the CG, and you can decide how you want to bin your regions:
freq = dinucleotideFrequency(SEQ,as.prob=TRUE)
hist(freq[,"CG"],br=100)


Simulating DNA sequence evolution in R I hope someone can lend their thoughts on the below code to generate DNA sequences under the Kimura-2-Parameter model of DNA substitution.

The issue is that each time the code is run and the haplotype distribution is examined, there always is a very skewed distribution.

What I would like is for a different distribution to be generated each time the code is run. that is, there should be multiple specimens sharing the same haplotype (instead of most sharing the most common haplotype).

Here is my code:

library(pegas)
set.seed(17)

sim.seqs <- TRUE
length.seqs <- 500
num.seqs <- 100 # number of DNA sequences
subst.model <- "K80" # nucleotide substitution model
transi.rate <- 1e-4 # transition rate
transv.rate <- transi.rate / 2 # transversion rate

if (sim.seqs == TRUE) {

  nucl <- as.DNAbin(c('a','c','g','t'))
  res <- sample(nucl, size = length.seqs, replace = TRUE, prob = rep(0.25, 4))

 if (subst.model == "K80") {
    transi.set <- list('a' = as.DNAbin('g'),
                       'c' = as.DNAbin('t'),
                       'g' = as.DNAbin('a'),
                       't' = as.DNAbin('c'))

    transv.set <- list('a' = as.DNAbin(c('c', 't')),
                       'c' = as.DNAbin(c('a', 'g')),
                       'g' = as.DNAbin(c('c', 't')),
                       't' = as.DNAbin(c('a', 'g')))

    transi <- function(res) {
      unlist(transi.set[as.character(res)])
    }

    transv <- function(res) {
      sapply(transv.set[as.character(res)], sample, 1)
    }

    duplicate.seq <- function(res) {
      num.transi <- rbinom(n = 1, size = length.seqs, prob = transi.rate) # total number of transitions

      if (num.transi > 0) {
        idx <- sample(length.seqs, size = num.transi, replace = FALSE)
        res[idx] <- transi(res[idx])
      }

      num.transv <- rbinom(n = 1, size = length.seqs, prob = transv.rate) # total number of transversions

      if (num.transv > 0) {
        idx <- sample(length.seqs, size = num.transv, replace = FALSE)
        res[idx] <- transv(res[idx])
      }

      res
      }

    }

  res <- matrix(replicate(num.seqs, duplicate.seq(res)), byrow = TRUE, nrow = num.seqs)

  class(res) <- "DNAbin"

  # write.dna(res, file = "seqs.fas", format = "fasta")

  h <- sort(haplotype(res), decreasing = TRUE, what = "frequencies")
  rownames(h) <- 1:nrow(h)

}

Output

h

Haplotypes extracted from: res

Number of haplotypes: 5
     Sequence length: 500

Haplotype labels and frequencies:

 1  2  3  4  5
96  1  1  1  1

If the code is run multiple times (without the seed), the same pattern emerges: h is always skewed toward the most dominant haplotype.

For example, I want to be able to run the code (without seed) and get something like

h

Haplotypes extracted from: res

Number of haplotypes: 5
     Sequence length: 500

Haplotype labels and frequencies:

  1    2    3    4   5
  35  25   20   15   5

i think I would have to specify a distribution for the haplotypes, but the number of haplotypes generated by the mutation process is not known beforehand.

Any thoughts?

What are you expecting these lines to do?

num.transi <- rbinom(n = 1, size = length.seqs, prob = transi.rate)
...
num.transv <- rbinom(n = 1, size = length.seqs, prob = transv.rate)

Using the variables that are specified in your code, they are extremely zero-skewed. This is not too surprising given the mutation rate and sequence length.

I expect that increasing the sequence length, or increasing the transition rate, should result in more haplotypes being generated.

After running this through myself, I see that increasing the mutation rate does increase the number of haplotypes, but the skew is retained. This is because each mutation creates a new haplotype, following something like an infinite allele model.

I guess this is not particularly surprising. It's unlikely that two random mutations from a particular haplotype will result in the same haplotype, and therefore any mutations (from any existing haplotype) will generally create something new.


How to download FASTA sequences from NCBI using the terminal? I have following accession numbers of the 10 chromosomes of Theobroma cacao genome.

NC_030850.1
NC_030851.1
NC_030852.1
NC_030853.1
NC_030854.1
NC_030855.1
NC_030856.1
NC_030857.1
NC_030858.1
NC_030859.1

I need to download these FASTA files using the terminal because I'm working on a server. I'll be grateful if someone can help me.
Thank you.

Alternatively, you can use the NCBI Entrez Direct UNIX E-utilities

Basically, you have to download the install file here: https://www.ncbi.nlm.nih.gov/books/NBK179288/bin/install-edirect.sh

In the terminal, install it using: source ./install-edirect.sh

Then, you can download your sequence by doing:

esearch -db nucleotide -query "NC_030850.1" | efetch -format fasta > NC_030850.1.fasta

And you should find your fasta sequence downloaded.

As you have several sequences to download, I think it will be quite easy to add this command into a little bash script to process all of them.

More readings about that: http://bioinformatics.cvr.ac.uk/blog/ncbi-entrez-direct-unix-e-utilities/ and https://www.ncbi.nlm.nih.gov/books/NBK179288/


Why did expression based subtypng of breast cancer gain much more acceptance than others This is may not be entirely technical question but rather a academic question. But the technique behind the application is within the scope of bioinformatics. So I would try to ask here that:
In each cancer type, there have been tons of papers that describe RNA-Seq based subtyping resulting in molecularly distinct subgroups with survival difference. However, only very few of them gained acceptance at clinical setting, such as breast cancer (luminal subtype A, luminal subtype B, normal breast-like subtype, HER-2 overexpression subtype and basal-like subtype).
My question is: why did this molecular subtyping of breast cancer gain much acceptance but not others? Did this one own some advantages that other did not?

In two words: incidence and funding
I'm not an expert on this topic, but I assume it has something to do with the incidence of breast cancer itself:

Breast cancer is the most common cancer in American women, except for skin cancers. Currently, the average risk of a woman in the United States developing breast cancer sometime in her life is about 13%. This means there is a 1 in 8 chance she will develop breast cancer. This also means there is a 7 in 8 chance she will never have the disease.

This incidence most likely amplifies how big of a problem society/funders/governments consider breast cancer to be. For example, a study performed by Northwestern University’s Feinberg School of Medicine and Robert H. Lurie Comprehensive Cancer Center showed that breast cancer research is by and large one of the most well-funded areas of cancer research.

This higher amount of funding leads to more and better resources and research into every aspect of breast cancer, from identification to treatment. I would wager better funding has had the effect of increasing the level of research performed in the usage of sub-typing for identifying people at risk of getting breast cancer. Given positive results and lots of literature backing up such a method as useful, its levels of acceptance in clinical settings increased.


Protein ligand docking: how to convert .pdb to .maps.fld? Hello I'm helping to develop a cloud docking tool for screening compounds, similar to Swissdock but with mass throughput and GPU optimizations.
Specifically helping screen existing drugs against coronavirus proteins (couple structures solved, rest TASSER).
I'm planning to use GPU optimized AutoDock-GPU , which takes in .maps.fld

I know you can use autogrid to select the bounding box and generate the .maps.fld, but I've been unable to figure out the workflow.
for preliminary screening I want to search the whole protein without specifying a bounding box.
Is there a script for converting protein.pdb to .maps.fld?

Collaborators welcome.
Thanks!

There are a few steps in between pdb and .maps.fld.
Here is the list of several scripts that can do tasks for you that you downloaded with autodock MGLTools: http://autodock.scripps.edu/faqs-help/faq/where-can-i-find-the-python-scripts-for-preparing-and-analysing-autodock-dockings. Look at the prepare_ files. Also note that the scripts come with a pre-bundled python, but you can install these scripts with conda.
Here are the basic order of things:
prepare_receptor4.py -r protein.pdb
prepare_ligand4.py -l ligand.mol2
prepare_gpf4.py -l ligand.pdbqt -r protein.pdbqt -y
autogrid4 -p protein.gpf
prepare_dpf4.py -l ligand.pdbqt -r protein.pdbqt

partial charge
Your PDB needs partial charges so first you convert it to PDBQT file.
This is true for your protein and your ligand. That is no complicated parameterisation for your ligand —which is not necessarily a good thing!
In this step it is essential to correct protonation if absent each have their own flags.
grid
The you make a grid file (gpf), which contains your grid. To make a box that spans the whole protein simply use the unit cell dimensions form the PDB. For an mmCIF, if you open the header dictionary you need _pdbx_struct_oper_list.matrix using MMCIF2Dict from biopython's PDB submodule. For a PDB file, just search for the CRYST1 line (cf. format).
Once you have have the grid you can run autogrid. This step will create the map.fld
Why a big box is bad for a VS
However, all those who take this approach opt for a box that is huge: as there is no solvent having a huge bounding has little penalty. However, this is considered a very bad strategy for a virtual screen. This not because it uses up computer resources as these are dirty cheap. Structural knowledge is important and targeting the active site of an enzyme stops it, while binding to the surface does nothing unless it's an interface. So manually reading what the stuff are and choosing the box wisely will save you a lot of time downstream in analysis. There is after all the saying "one week in the lab will save you an hour in the library"...
Protein and ligand
Lastly, merge the protein with the small molecule in a dpf file.
Caveat against docking with models
It is not a good idea dock against models. Swissmodel (unless virtually identical), I-Tasser, Phyre or EVFold models. Docking is very sensitive to small structural difference that may have resulted from the modelling, so docking to models is highly discouraged. For Coronavirus protein, dock against SARS, which is highly similar. Or at the very least pay strong attention to the I-Tasser C-scores and discard the bottom 2/3 of their models.
Note about coronavirus
Your objective of docking coronavirus may be a bit too late however.
I did a wee interactive summary of the literature about the solved structure which is the protease. Proteases are very easy to rationally design drugs for.
There are already:

The Covid Moonshot project is an open collaboration between a fragment-screening X-ray facility (Diamond XChem) and many researchers from various disciplines (myself included).
a paper that has an empirically validated coronavirus specific ligand
many many manuscripts about virtual screens in bioarxiv/chemarxiv.
licensed HIV protease inhibitors lopinavir and ritonavir that were initially suspected to be effective —a clinical trial

Also you may want to check out Galaxy project, a tool to run pipelines for genomics, but have started moving into biochemistry from genomics. Eg. https://covid19.galaxyproject.org/cheminformatics/.


Where did the BIND website go? I was trying to reach the BIND (Biomolecular Interaction Network Database) but the website is down.

Where can I find something similar? Has this platform moved somewhere or what?

Answer from @heathobrien converted from comment:
looks like that link has been dead for years:
https://www.biostars.org/p/16763/
Answer from @jason-aller converted from comment:
The biostars article mentions this:
http://baderlab.org/BINDTranslation


Download data from the Human Microbiome Project via ascp I have asked this question in biostars, but I am trying here as well as people working with "omics" data might be able to help. I think my issue relates understanding how large data storage on online server works. I am trying to download large data from the Human Microbiome Project. These are available on online server and can be downloaded via aspera (ascp).

I can't figure out the right path to get the data. On their website, they give an example of command to retreive an example dataset:

 ascp -l 100M -QT user@aspera.hmpdacc.org:/ReferenceGenomes/test_genome.tar /path/to/local/directory

I executed this, using my own username, from this I am asked my password, and I can authenticate. But then I get the following error:

(Error: Server aborted session: No such file or directory)

Then, I downloaded the manifest file via the HMP portal. In this manifest, they provide the path to which these files are supposed to be available. I tried to execute again the command with one of those path, provided by the manifest file:

ascp -l 100M -QT user@aspera.hmpdacc.org:/files/5a950f27980b5d93e4c16da124000378  /path/to/local/directory

I again get the same error of "No such file or directory".

The full url provided by the manifest

I tried putting different chunks of this address into to ascp command but nothing works. I think there is something I don't understand about how to access the server via aspera and what is expected to go in this path. When I copy paste this url in google, this directly download the file I want. So I'm close, but I need to do it via command line as I need to download hundreds of these files and they are huge (6-10GB each).

How can I download the data of the human microbiome project using ascp?

I can't test for ascp because they do not longer support my current setup but you can follow the FTP instructions. With FileZilla:

url: public-ftp.ihmpdcc.org