Extracting genes differentially expressed by Wilcox test I have 2 cell types (NOF and CAF) cultures individually and also have been co-cultured with tumour like this picture.

https://ibb.co/x8Ty1Wz

If this is the log cpm of these 4 samples

> head(onco[,1:4])
             A10       A11        A12        A2
A2M    11.755100 12.145143 12.9787977 12.407221
ABCB11  2.955457  2.278363  0.1229266  1.493507
ABCC2   6.711643  6.974025  0.1229266  7.628936
ABCC6   6.443439  6.094670  6.1389496  7.483394
ABCF1   9.995098 10.045613  9.6627424 10.314842
ABCG2   7.125483  7.448127  6.8669702  8.492895
>

I have used log of count per million (cpm) for Kruskal-Wallis rank sum test that says these samples are significantly difference like below;

> NOF1=c(NOF1$`onco$A10`)
> NOF2=c(NOF2$`onco$A11`)
> CAF1=c(CAF1$`onco$A12`)
> CAF2=c(CAF2$`onco$A2`)
> my_data <- data.frame(
+     group = rep(c("NOF1", "NOF2","CAF1","CAF2"), each = 719),
+     cell_type = c(NOF1,NOF2,  CAF1,CAF2)
+ )

> head((my_data))
  group cell_type
1  NOF1 11.755100
2  NOF1  2.955457
3  NOF1  6.711643
4  NOF1  6.443439
5  NOF1  9.995098
6  NOF1  7.125483
>

> kruskal.test(cell_type ~ group, data = my_data)

    Kruskal-Wallis rank sum test

data:  cell_type by group
Kruskal-Wallis chi-squared = 22.006, df = 3, p-value = 6.505e-05

>

I then did Pairwise comparisons using Wilcoxon rank sum test to see which samples are different with each other.

> pairwise.wilcox.test(my_data$cell_type, my_data$group,
+                      p.adjust.method = "BH")

    Pairwise comparisons using Wilcoxon rank sum test

data:  my_data$cell_type and my_data$group

     CAF1   CAF2   NOF1
CAF2 0.0058 -      -
NOF1 0.2809 4e-05  -
NOF2 0.6818 0.0016 0.3811

P value adjustment method: BH
>

My Question is; how I can extract the genes behind the difference between each pair of samples? I mean something like differential expression by non parametric test. also I am not sure if I am right in using log cpm values for these test

If possible give me a hand to correct myself and get intuition

Thanks

All of the genes are behind the difference, you're doing a rank test comparing all of them against each other at once. Your results have more to do with samples having slightly different outlier genes than any interesting biology. You don't have enough replicates to do a non-parametric test.


Understanding SNP coding for association analysis I'm working on a project about detecting SNP association with a disease. As I understand, SNP is a single variation of the nucleotide that occurs for more than 1% of the population. So, if a gene is ATTG and there's a variation ATAG then the SNP is A. However, I couldn't connect this idea with the dataset in hand. The rows in my dataset represent each patient and the columns contain SNP information. For example:

ID  exm355 exm615
1     T_T   A_C
2     T_T   C_C
3     A_T   C_C

I have no idea why the SNP columns contain 2 nucleotides(T_T, A_T, A_C, C_C). As the definition of SNP, I thought it should show only the variant nucleotide or am I misunderstanding anything? How could I interpret T_T or C_C and how could I know which nucleotide is the variation from the common ones in the population?

The source says that the case and control subjects were genotyped by the Australian Genome Research Facility using the Illumina HumanExome BeadChip-12v1_A, which exonic content consists of $>$ 250 000
markers representing diverse populations and a range of common
conditions. An individual was duplicated to test genotyping reliability and
quality. The identity by descent (IBD) matrix between all pairs of individuals
was estimated after linkage disequilibrium (LD) pruning and used for
quality control and for the mixed linear models analyses.

Humans are diploid. This means they have two copies of each locus (with the exception of sex chromosomes in males). Your two alleles for each individual represent the two copies of this variant locus.


HISAT2 - are all multi mapped reads reported by default? While going through the HISAT2 manual, one can see this:

  Two alignments for the same individual read are "distinct" if they map
  the same read to different places. Specifically, we say that two
  alignments are distinct if there are no alignment positions where a
  particular read offset is aligned opposite a particular reference
  offset in both alignments with the same orientation. E.g. if the first
  alignment is in the forward orientation and aligns the read character
  at read offset 10 to the reference character at chromosome 3, offset
  3,445,245, and the second alignment is also in the forward orientation
  and also aligns the read character at read offset 10 to the reference
  character at chromosome 3, offset 3,445,245, they are not distinct
  alignments.

  Two alignments for the same pair are distinct if either the mate 1s in
  the two paired-end alignments are distinct or the mate 2s in the two
  alignments are distinct or both.

  Default mode: search for one or more alignments, report each

and then there is -k parameter (analogous to bowtie2 -K parameter)

which defines:

  HISAT2 searches for up to N distinct, primary alignments for each
  read, where N equals the integer specified with the -k parameter.
  Primary alignments mean alignments whose alignment score is equal or
  higher than any other alignments. It is possible that multiple
  distinct alignments have the same score. That is, if -k 2 is
  specified, HISAT2 will search for at most 2 distinct alignments. The
  alignment score for a paired-end alignment equals the sum of the
  alignment scores of the individual mates. Each reported read or pair
  alignment beyond the first has the SAM 'secondary' bit (which equals
  256) set in its FLAGS field. See the SAM specification for details.

This is a bit confusing, as the first quote means that HISAT2 finds all distinct valid alignment and reports them by default, and the second quote says that when -k option is used, it will only report a subset of alignments which it would have reported by default. Am I correct in this inference or does one needs to use -k option to enable "multi-mapping"?

You said:

  ...the first quote means that HISAT2 finds all distinct valid alignment and reports them by default.

Not exactly. The first section you quoted doesn't say it will find and report all distinct alignments, but rather that it would find and report one or more distinct alignments. It's a subtle difference, but an important one.

That section describes how distinct mappings are defined, not how many are reported. The second section you quoted, as well as the documentation for the -k parameter, are all consistent with this interpretation.

-k <int>        It searches for at most <int> distinct, primary alignments for
                each read. Primary alignments mean alignments whose alignment
                score is equal or higher than any other alignments. The search
                terminates when it can't find more distinct valid alignments, or
                when it finds <int>, whichever happens first. The alignment
                score for a paired-end alignment equals the sum of the alignment
                scores of the individual mates. Each reported read or pair
                alignment beyond the first has the SAM 'secondary' bit (which
                equals 256) set in its FLAGS field. For reads that have more
                than <int> distinct, valid alignments, hisat2 does not guarantee
                that the <int> alignments reported are the best possible in
                terms of alignment score. Default: 5 (HFM) or 10 (HGFM)

                Note: HISAT2 is not designed with large values for -k in mind,
                and when aligning reads to long, repetitive genomes large -k can
                be very, very slow.

By default, 5 or 10 distinct mappings are reported (depending on the mode), but you can control this with the -k parameter.


Problem: "Pair-end" reads scRNA seq data (Drop-seq) In case of Drop-seq, we have paired end data.

Read 1: Cell code + UMI (unique molecule identifier)

Read 2: The transcript information

But I have a problem/doubt with the sample I am working on.

The sample I am using is the following (Check the "Reads" tab).

As you know the Drop-seq is "paired-end", we are expected to see two reads for a spot. Although this sample say paired-end, it has only one read per spot.

For example I can share a link of a different scRNA-seq data where you can properly see two reads for a spot

Example sample, you need to check the "Reads" tab.

Where I am going wrong?

I asked one of the main authors of the paper. The following is the reply I got :

  I recommend that you download the aligned BAM files that are hosted in
  the same GEO record. Read 1 is already processed into the cell and UMI
  barcodes and held as custom tags (XC and XM) in the BAM files. The
  cells are already barcode-corrected, so if you use those files, your
  cell barcodes will line up with mine; if you start from FASTQs, they
  will not. For most aligners, you can just use the BAM file as input to
  realign. (It has all reads, even those that did not align.)

But I could not find any "XM" or "XC" keywords in the bam file :(

To understand his reply you have to be familiar the processing steps of the Drop-seq.

Looks like they have submitted some kind of processed data. I could not figure out how much the data is processed. I am trying to use the data starting at different processing steps but I could not figure out how much the data is processed.

I received the following message from the author of another Drop-seq paper:

  After much conversation with the folks at the SRA, it appears the SRA
  compression algorithm automatically deletes certain subsets of bam
  tags and does so in cryptic way (in the sense of not alerting the
  uploader). Obviously this is going to be an issue for many single cell
  datasets encoded this way, and I hope the SRA will eventually change
  their system to accommodate.

If you look at the "Download" tab of the same page, there now exists a link to download the originally uploaded BAM file.

Alternately, you can use the SRA AWS tutorial to download the files on that platform.


Making a plot cleaner By this code I produced this tSNE

> length(expr)
[1]  1 209
> aaa=as.matrix(signif(expr,digits=1))
> colnames(aaa)="value"
> aaa=as.data.frame(aaa)
> cc <- AddMetaData(object = seurat, metadata = aaa)
> cc <- SetAllIdent(object = cc, id = "value")
> head(as.matrix(cc@meta.data))
     nGene  nUMI      orig.ident      res.0.7 CELL   STAGE GENO dataset
s1.1 "4331" " 373762" "SeuratProject" "0"     "s1.1" "H16" "WT" "1"
s1.2 "5603" "1074639" "SeuratProject" "0"     "s1.2" "H16" "WT" "1"
s1.3 "2064" "  49544" "SeuratProject" "0"     "s1.3" "H16" "WT" "1"
s1.4 "4680" " 772399" "SeuratProject" "1"     "s1.4" "H16" "WT" "1"
s1.5 "3876" " 272356" "SeuratProject" "1"     "s1.5" "H16" "WT" "1"
s1.6 "2557" " 122314" "SeuratProject" "0"     "s1.6" "H16" "WT" "1"
     stage.nice celltype value
s1.1 "H16"      "0"      "-2e-01"
s1.2 "H16"      "0"      " 3e-02"
s1.3 "H16"      "0"      "-1e-01"
s1.4 "H16"      "1"      "-4e-02"
s1.5 "H16"      "1"      " 2e-01"
s1.6 "H16"      "0"      "-8e-02"
> TSNEPlot(object = cc)

As you are seeing the range of expression still is very sporadic or when simply use seurat featueplot, again plot is not informative

> feature.plot(seurat,features.plot ="DDB_G0267412")

I fed my scaled data from seurat to single cell orange software, when I plotted one gene on tSNE gave me this

You are considering how much the plot is clean

How can I  have such a clean plot please? I mean how I can obtain such a range in colouring my cells

> head(cc@meta.data)
     nGene    nUMI    orig.ident res.0.7 CELL STAGE GENO dataset stage.nice
s1.1  4331  373762 SeuratProject       0 s1.1   H16   WT       1        H16
s1.2  5603 1074639 SeuratProject       0 s1.2   H16   WT       1        H16
s1.3  2064   49544 SeuratProject       0 s1.3   H16   WT       1        H16
s1.4  4680  772399 SeuratProject       1 s1.4   H16   WT       1        H16
s1.5  3876  272356 SeuratProject       1 s1.5   H16   WT       1        H16
s1.6  2557  122314 SeuratProject       0 s1.6   H16   WT       1        H16
     celltype  value
s1.1        0  0.080
s1.2        0  0.100
s1.3        0 -0.060
s1.4        1  0.007
s1.5        1 -0.060
s1.6        0  0.010
>

It doesn't look like this data has been clustered, and you are plotting the 'value', rather than the cluster. If you want to represent values on a tSNE using Seurat, a better function to use is the FeaturePlot, which will colour the points in a more meaningful colour gradient based on their value.
[It's possible this function wasn't available at the time this question was written]
If, instead, you want to plot clusters, then you should carry out some clustering on the data first using FindNeighbors/FindClusters.
There's a great introductory tutorial to guide you into using Seurat on the Satija lab page:
https://satijalab.org/seurat/articles/pbmc3k_tutorial.html


Are mitochondrial genes to exclude in scRNA-seq such as ribosomal genes? In this answer, it is stated that ribosomal genes should be excluded prior to normalization in scRNA-seq as contaminants.

Do mitochondrial genes have to be excluded as well? I plotted the top 50 expressed genes for a specific dataset and they tend to appear often (for example MT-ATP6). My assumption is that, given that they work for mitochondrial function and may be highly expressed, they can dilute the signal of gene differential across cell types but expressed at lower levels. Is this biologically sound?

Additionally, in this course, mitochondrial genes are used to filter cells when they contribute above a certain ratio of total RNA of single cells. However, I could not find an explanation for this filtering. What is the underlying process this procedure is used to account for?

According to Ilicic et al. (2016), on upregulation of mtRNA in broken cells:

  There is an extensive literature on the relationship between mtDNA, mitochondrially localized proteins, and cell death [34, 35]. However, upregulation of RNA levels of mtDNA in broken cells suggests losses in cytoplasmic content. In a situation where cell membrane is broken, cytoplasmic RNA will be lost, but RNAs enclosed in the mitochondria will be retained [...]

According to this, I think that mtRNA levels relative to endogenous RNA can be used as a control for low quality (broken) cells, in a similar way to ERCC spike-ins.

Ilicic T, Kim JK, Kolodziejczyk AA, Bagger FO, McCarthy DJ, Marioni JC, Teichmann SA. Classification of low quality cells from single-cell RNA-seq data. Genome biology. 2016 Feb 17;17(1):29.


Gene set enrichment analysis on differential phosphorylation sites I have:

A list of differentially phosphorylated sites in a knockout condition. Some genes contain as many as 70 possible phosphorylation sites; others contain only one.
A list of genes belonging to a specific gene set annotation.

How can I test the differentially phosphorylated proteins for enrichment of this annotation?

A few ideas I’ve considered:

Ignore the number of phosphorylation events detected within a gene and simply count the gene as differentially phosphorylated if it contains at least a single site that is differentially phosphorylated. Compare the enrichment score for this selected set to the enrichment score for the set genes containing at least one site that is not differentially phosphorylated. The problem here is that genes with a large number of phosphorylation sites have almost no influence in the enrichment score, since it’s almost certain that they’ll have at least one site that is differentially phosphorylated and at least one that is not.
Mark each candidate phosphorylation site as either “in the set” or “not in the set” based on the protein in which it’s found. Then perform the enrichment analysis using the set of annotated phosphorylation sites instead of using the traditional enrichment analysis performed at the gene level. The potential problem with this approach is that it may place too much influence on genes with many potential phosphorylation sites.
Aggregate all candidate phosphorylation sites within a gene and use some numerical threshold to determine whether the gene is differentially phosphorylated or not. (There are various ways that this could be done.) Then perform enrichment analysis using the resulting set of differentially phosphorylated genes. A possible problem here is that some of the phosphorylation sites may be more functionally important than others, so it’s not clear how to weight the relative importance of individual phosphorylation sites within a gene.

I realize the goal here isn’t well-defined mathematically; I’m mainly curious what approach makes the most sense given the biological context. Of the above approaches, I’m currently leaning towards approach 2 because it’s straightforward to implement and it at least attempts to account for the variable number of phosphorylation events within a gene.

NOTE: I also have normalized phosphoprotein and protein abundances for all of these sites, obtained from mass spec. So if the solution requires an alternative method of computing differential phosphorylation, that’s fine.

Your dataset looks perfect for the SUMSTAT enrichment test.

You need to come up with a statistic representing your gene. The simplest ideas here are the number of sites or, probably better, a proportion of phosphorylated sites.
Now you can compute a statistic for every gene set, for example just sum (SUMSTAT). You can also have average, sum of squares or something else.
Geta a null distribution of your gene set statistics by permutations. You can either randomize statistics for genes, or just randomly assign phosphorylated sites across the genome keeping the total number of sites.

Now you can compute p-value just by comparing your value with the null-distribution.

You have to be cautious of two things:

You are testing multiple gene sets, i.e. you perform a statistical test for every gene set. Therefore correct your significance for multiple testing. I suggest using FDR for this.
Beware of potential biases when doing permutations. The obvious one is gene length, i.e. longer genes have higher chance of getting at least one site. But there can be other ones like GC-content, chromosomes, etc. You can overcome this by using more realistic permutations, or by controlling for the potential correlations. You can get a few ideas on fighting biases from this paper (sorry for the self-advertisement).


Seurat FindMarkers() output, percentage I have generated a list of canonical markers for cluster 0 using the following command:

cluster0_canonical <- FindMarkers(project, ident.1=0,  ident.2=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14), grouping.var = "status",
min.pct = 0.25, print.bar = FALSE)

And I got the following output, but I do not understand why the canonically expressed ribosomal markers seem to be expressed at high levels in both the cells within the cluster AND cells that are NOT in the cluster (i.e. both pct.1 and pct.2 are close to 1, but P value is highly significant). We would not expect this for canonical markers - we would instead expect the genes to be expressed at high levels in the cells within the cluster but at low levels in cells not within the cluster. I don't seem to have such an issue with non-ribosomal genes within cluster 0 and cluster 0 is the only cluster where I get a high number of ribosomal genes as canonical markers.

Your interpretation of pct.1 and pct.2 is misguided. If pct.1 == 1, the given gene was detected in all of the cells. This carries no information about high expression or low expression in the cell, just that it was detected. avg_logFC is the indication of expression levels. So your expectation, that

  "we would instead expect the genes to be expressed at high levels in the cells within the cluster but at low levels in cells not within the cluster"

holds true, as indicated by the positive avg_logFC and significant p-values of the listed RP genes.

Seurat use nautral log, so the FC of RPS6 in cluster 0 vs. all other clusters indicated is 2.718281828459^.55947=1.750.

Seurat FindMarkers() documentation.


For a data.frame named foo:

sink('some_output_file.txt')
uvals = unique(foo[,1])
for(uIDX in c(1:length(uvals))) {
    gname = uvals[uIDX]
    IDX = which(foo[,1] == gname)
    cat(sprintf("%s\t%s\n", gname, paste0(foo[IDX,2], collapse="\t")))
}
sink()


differential expression multiple studies I was wondering if we can perform differential expression on samples from different library and read lengths.

What we intend to do is,

n = 6 tumors (recently sequenced using Hiseq 2500 read length 150, paired end)
n = 5 normals (previously sequenced using Hiseq read length 100, paired end, these are not matched normal)
n = 100 ( normals, from tcga read length 100, paired)

What we would like to do is compare n=6 Vs 100 + 5 normals. What steps should we be considering before we move ahead with this analysis?

Adding experiments 1 & 2 should be fine without modification on the computer side of things, as long as your sample preparation is consistent. To be ultra-cautious, you could add in the experiment as an additional covariate for the differential expression model.

I'm less confident about adding in TCGA data; it's probably going to look better than your own samples, which might cause problems in calculating differential expression from raw reads. You can do a separate normalisation/correction of the TCGA data, but that will limit the inferences that can be made about differential expression -- without raw count information, the shot noise associated with read sampling cannot be included in the statistical model. I notice that TCGA has already done their own filtering / expression, and that would be preferable to use (mRNA sequencing, data subtype 'Gene'):

https://tcga-data.nci.nih.gov/docs/publications/tcga/datatype.html

However, even with data that has been properly normalised and filtered to exclude batch effects, there will still be population differences that influence results. My recommendation would be to do some subsampling to reduce the population effects, something like at least the following comparisons:

6 Vs 100 + 5 normals
6 Vs 100 normals
6 Vs 50 + 3 normals [randomly sampled from the two normal groups]

And only consider differential expression that is consistently observed in all three of these comparisons. Note that the third "subsampled" dataset can be generated more than once, taking different subgroups each time.

If you are going to use separately normalised data, then DESeq2 may not be appropiate (because it expects count-level data). Something like Limma could work for DE analysis, see chapter 15 of the user's guide for RNASeq analysis. Unfortunately I'm not very familiar with using Limma, so can't provide any advice on things like what to watch out for in the input data and/or results.


What are the de facto required fields in a SAM/BAM read group? The SAM specification indicates that each read group must have a unique ID field, but does not mark any other field as required.

I have also discovered that htsjdk throws exceptions if the sample (SM) field is empty, though there is no indication in the specification that this is required.

Are there other read group fields that I should expect to be required by common tools?

The sample tag (i.e. SM) was a mandatory tag in the initial SAM spec (see the .pages file; you need a mac to open it). When transitioned to Latex, this requirement was mysteriously dropped. Picard is conforming to the initial spec. Anyway, the sample tag is important to quite a few tools. I would encourage you to add it.


Why is this makeblastdb command not working? I am trying to make local databases for the ncbi-blast+ package (version 2.60). I am doing so for 4 T-cell receptor genes. 3 of the 4 (TRAV, TRAJ, TRBV) have worked fine, but I am having problems with TRBJ. I say this because they should be the same format and if one works the rest works.

They take the format below:

>K02545|TRBJ1-1*01|Homo sapiens|F|J-REGION|749..796|48 nt|3| | | |15 AA|15+0=15| | |
NTEAFFGQGTRLTVV

I used the following command (like I say, it was fine for the rest, oddly enough):

makeblastdb -in TRBJ_AA.fa -parse_seqids -dbtype prot

and I get the following error message:

BLAST options error: TRBJ_AA.fa does not match input format type, default input type is FASTA

One issue I checked was the length of the sequences being too short, but I believe 11 is the minimum length required (for inputs, don't know if that is the same for databases)- but would someone be able to confirm, that strikes me as the most likely candidate?

I'm happy to use some other metric, but I was hoping on keeping BLAST for consistency's sake for the time being.

Thank you.

edit: I know it's to late, but for reference I'm on Ubuntu 16.04.2 LTS x86_64 for all steps.

This looks like a bug in makeblastdb. Removing the final | | | from your sequence's description makes it work:

>K02545|TRBJ1-1*01|Homo sapiens|F|J-REGION|749..796|48 nt|3| | | |15 AA|15+0=15

So does removing everything after the 1st space:

>K02545|TRBJ1-1*01|Homo sapiens|F|J-REGION|749..796|48

There is nothing in the FASTA definition that would explain this and, even more weirdly, the presence of the | | | alone isn't enough to trigger the issue. I played around with this a little and it looks like it doesn't like the combination of + or -along with empty | |, but there's no clear pattern I can see:

These fail:

>| | |15-| | |
>| | |15+| | |
>| | |15-1| | |
>| | |15+1| | |

These work:

>| |15+| |
>| | |151| | |

So yeah, sounds like a bug. I suggest you report it to the NCBI devs. In the meantime, a simple fix would be to just keep the first space-separated field of your fasta header. Pass your fasta file through this:

awk '{print $1}' fasta.fa > new.fa

That will print the 1st field of each line and, assuming you don't have spaces in your sequence, it should effectively sort this out for you.


What is "block-compressed" file in samtools? SAMtools returned an error message for my gzipped genome FASTA:

  Indexed block-compressed FASTA file cannot be handled

The source code for the message is here.

// check if the it is a valid block-compressed file
if (IOUtil.isBlockCompressed(path, true)) {
    throw new SAMException("Indexed block-compressed FASTA file cannot be handled: " + path);
}

What's a "block-compressed" file?

A block compression usually refers to compressing your file into a series of small blocks (with a tool like bgzip). This allows indexing in that the index can record which record lives in which block so that the whole file does not need to be decompressed and therefore can be accessed randomly. This isn't true of files compressed as a single block using gzip.

As far as I was aware samtools only works with uncompressed FASTA files.


move files of a different group to separate directory EDITED

I have list of my samples in first column and corresponding file name for my .vcf files in second column of a file name clin_name.txt like below. For example name of one vcf is LP6008031-DNA_A04_vs_LP6008032-DNA_D03.snp.pass.vcf

OC/AH/120/E/T3  LP6005409-DNA_F01
OC/AH/126/E/T1  LP2000325-DNA_A01
OC/AH/136/E/T3  LP6005409-DNA_E02

This is the link of my file https://www.dropbox.com/s/mvzf47feiqe7x69/clin_name.txt.txt?dl=0

In the second file name pre_post.list I have information of pre or post treatment for each sample like this

pre-treated post-treated
OC/AH/120   OC/AH/065
OC/AH/126   OC/AH/069
OC/AH/136   OC/AH/074
OC/AH/160   OC/AH/084

This is link of this https://www.dropbox.com/s/nhv3odjkdjdp7hz/pre_post.list.txt?dl=0

For moving different .vcf files based on the pre or post treatment to separate folders I typed

cat pre_post.list | awk '{print $2}' | sed 's@.*/@@'  > post_list.txt
cat pre_post.list | awk '{print $1}' | sed 's@.*/@@'  > pre_list.txt

Then I have

pre-treated
120
126
136
160
182

post-treated
065
069
074
084

I then typed ./test.sh in which I wrote

for item in  $(cat  pre_list.txt | grep -Eo '[0-9]{1,4}');
    do if grep -q $item clin_name.txt ;
    then out=$(grep $item clin_name.txt | awk  '{print $1}');
         cp *$out* /temp/hgig/fi1d18/TRG45/snp/
    else echo post_list $item;  fi;  done;

That returns

cp: cannot stat `*OC/AH/069/S/T3*': No such file or directory
[fi1d18@cyan01 snp]$

Or I run

for item in  $(cat  post_list.txt | grep -Eo '[0-9]{1,4}');
    do if grep -q $item clin_name.txt ;
    then out=$(grep $item clin_name.txt | awk  '{print $1}');
         cp *$out* /temp/hgig/fi1d18/TRG45/snp/pre/
    else cp *$out* /temp/hgig/fi1d18/TRG45/snp/post/;  fi;  done;

That returns

cp: cannot stat `OC/ST/026/S/T1*': No such file or directory
[fi1d18@cyan01 snp]$

I really don't know how to solve this :(

You need to bring together the information in the two tables, I prefer R for such actions. The common columns between the two tables do not quite contain the same information hence additional data cleaning was required. Moreover the table with treatment information had to be changed to long format with melt. And I have filtered the merged_data so as to conation only "pretreatment" info but you can add a if/else statement to the for loop to control how/where you will move your files.

library(data.table)

clin_name <- fread("clin_name.txt", header = TRUE)
clin_name$sample_id <- sub("/E/T.", "", clin_name$sample_id)

pre_post <- fread("pre_post.list", header = TRUE)
print(pre_post)
pre_post <- melt(pre_post,
     measure.vars = c("pre-treated", "post-treated"),
     variable.name = "treatment",
     value.name = "sample_id")

> print(pre_post)
      treatment sample_id
1:  pre-treated OC/AH/120
2:  pre-treated OC/AH/126
3:  pre-treated OC/AH/136
4:  pre-treated OC/AH/160
5: post-treated OC/AH/065
6: post-treated OC/AH/069
7: post-treated OC/AH/074
8: post-treated OC/AH/084

merged_data <- merge(clin_name, pre_post)
print(merged_data)

pre <- merged_data[merged_data$treatment == "pre-treated", "file_name"]

> print(merged_data)
   sample_id         file_name   treatment
1: OC/AH/120 LP6005409-DNA_F01 pre-treated
2: OC/AH/126 LP2000325-DNA_A01 pre-treated
3: OC/AH/136 LP6005409-DNA_E02 pre-treated

# invoking a terminal command to move each file to a directory called "pre"
for (i in pre$file_name) {
  command <- paste0("mv ", i, " pre/")
  system(command)
}


Genome databases of microbial eukaryotes I am interested in collecting the reference genomes of marine protists and would like to know which databases can be used for such tasks, and is there a single, centralized database of protistological genomic data?

You can find many protist genomes on ensembl

http://protists.ensembl.org/index.html


Weighted Network Analysis Book by Steve Horvath, taking the power of a matrix In Horvath's Weighted network analysis, the author mentions that M^b,  where b is a real number and M is a matrix, denotes the element-wise power. This seems odd to me. The way I was taught, M^b would mean multiplying the matrix M, b times. Should element-wise power not be M^(ob) where where o denotes the Hadamard power?

You CAN calculate M^b for b being float number, complex number or even another matrix.
As well as you can calculate exp(M), log(M), sin(M) whatever... That is very standard for math people.

Let me give two ways how you can think of it.

1) M^b = exp(b log(M)), so you need to define what is exp and log for matrices.
That can be done using power series: exp(M) = 1 + M + M^2/2! + M^3/3! + ...
log(1+M) = M - M^2/2 + M^3/3 - ...

2) More general way to calculate any function of matrix f(M) , looks like that:
if your matrix is diagonal, then it is true that: f(diag(d_i)) = diag( f(d_i)  ) .
So the question how to reduce any matrix to diagonal one.
That can be done for almost all matrices by conjugation M = C D C^{-1}.
Then f(M) = C f(D) C^{-1}.


sum of products of mapq and mapped bases for each read in a from a BAM file Given a BAM file, I'd like to calculate the sum, over all reads, of the mapping quality and the number of mapped bases (i.e. number of M's in the CIGAR string).

For example, given two reads like this:

read1: mapq=40, mapped_bases=10
read2: mapq=20, mapped_bases=50

The number I'm looking for is:

(40*10) + (20*50) = 1400

One option here is to stream the file and get the two quantities line-by-line. But this seems likely to be slow. Thoughts welcome.

BAM files are stored line by line, so your only option is to process it that way. In python with pysam:

import pysam

bam = pysam.AlignmentFile("something.bam")
tot = 0
for read in bam.fetch():
    if read.is_unmapped:
        continue
    tl = 0
    for op, l in read.cigartuples:
        if op in [0, 7, 8]:
            tl += l
    tot += tl * read.mapping_quality
bam.close()
print("the total is {}".format(tot))


Single-cell RNA sequencing (scRNA-seq): filtering cells by transcript counts, how to choose cutoffs? I am running a notebook with example for the MAGIC algorithm.

In the data preprocessing step, a filtering operation is required to filter out cells with a small count of transcripts.

My question is twofold:

Why filtering out cells with a small count of transcripts?
How to decide the cutoff, i.e. the minimum number of transcripts each cell should have to not be removed from the dataset?

I attach some plots from the original notebook, which display some distributions of the single-cell rnaseq dataset.

By reading this tutorial:

  Wells with few reads/molecules are likely to have been broken or failed to capture a cell, and should thus be removed.

This answers the first part of the question.

For the second part, i.e. how to decide cutoff, it still seems a bit arbitrary to me, but some authors trim the cells on the left tail of the highest mode, in order to remove other modes (likely due to noise) on the left tail, considering the distribution of cells dependent on the sum of reads (molecules per cell, first plot).


assigning different characters to present and absent genes from blast output using perl or python I have file called pangenome.txt which has only one column containing the names of the genes

head pangenome.txt
    Ha1_00044161
    Ha11_00031248
    Ha16_00021664
    Ha8_00042762
    Ha4_00025807

I have another bunch of files in a tabular blast format (330 in total: individual1.blast.txt,individual2.blast.txt, ...,individual330.blast.txt ) which are the outputs of Diamond protein aligner. Below is part of individual1.blast.txt file

head blast.individual1.txt
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha16_00021664   77.4    31  7   0   94  2   49  79  7.3e-159.3
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha8_00042762    74.2    31  8   0   94  2   1   31  9.6e-158.9
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha4_00025807    77.4    31  7   0   94  2   1   31  1.3e-058.5
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha16_00022939   77.4    31  7   0   94  2   295 325 2.1e-057.8
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha12_00032339   71.0    31  9   0   94  2   2   32  2.8e-057.4
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha16_00021553   71.0

I want to know which genes from the pangemome.txt is present in my individuals file and which are absent.
for the example above this is my desired output:

head output.individual1.txt
    Ha1_00044161   0
    Ha11_00031248  0
    Ha16_00021664  1
    Ha8_00042762   1
    Ha4_00025807   1

My desired out put is a file with 2 columns (tabular format), the first column the list of the all genes form pangemome.txt pan, and the second column is the presence and absence of these genes in the individuals files if the gene is present present in the blast.individual1.txt, I want to assign 1 and if it is absent I want to assign 0.
Can anyone help me to do this?
Thanks a lot

I wouldn't do this in Python, myself. This is a very simple text parsing problem and the standard *nix tools will be able to do it very easily. For example with awk:

$ awk 'NR==FNR{a[$1]=0; next}
       {
        if($2 in a){
            a[$2]=1
        }
       }END{
        for(i in a){
            print i,a[i]
        }
       }' pangenome.txt blast.individual1.txt
Ha1_00044161 0
Ha4_00025807 1
Ha11_00031248 0
Ha16_00021664 1
Ha8_00042762 1

The above is just a one-liner, I expanded it for clarity. You can run it on as many files as you like:

awk 'NR==FNR{a[$1]=0; next}{if($2 in a){a[$2]=1}}END{for(i in a){print i,a[i]}}' pangenome.txt blast.individual*.txt

For many large files, the fastest approach would probably be to just grep the names (unlike the awk, this will also match if the target gene is in the 1st field, but that's probably not an issue for you):

$ grep -wof pangenome.txt blast.individual1.txt.orig
Ha16_00021664
Ha8_00042762
Ha4_00025807

That will give you the list of genes found in the target blast files. To get the 1 and 0, you can do:

$ sort <(grep -owf pangenome.txt blast.individual1.txt) pangenome.txt |
    uniq -c | awk '{print $2,$1-1}'
Ha1_00044161 0
Ha11_00031248 0
Ha16_00021664 1
Ha4_00025807 1
Ha8_00042762 1

For multiple files, just run it on all of them:

$ sort <(grep -owf pangenome.txt blast.individual*.txt) pangenome.txt |
    uniq -c | awk '{print $2,$1-1}'


Database with antigens and current vaccines I am seeking a virological database comprising the following information:

vaccine.
antibody
receptor.

Preferably hosting multiple strains/samples of a given virus samples of strain.

Any assistance would be welcome.

Given "any virus will do", the viruses that are "head and shoulders" above the rest in terms of volumes of data are HIV and influenza.

The influenza database is here

The influenza viral spike that attaches to the cell receptor is the HA protein – hemagglutinin. The cell receptor is sialic acid – a small sugar that is attached to many different proteins on the cell surface. You can of course simply search PDB for HA or else Genbank.


Curl error when installing Bioconductor packages with Microsoft R Open 3.5.0 I'm trying to install Bioconductor packages with Microsoft R Open 3.5.0 seems to fail not sure what is the issue.

This is what I'm trying to do

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

But I keep getting an error about "curl" had non-zero exit.
How can I install the packages of Bioconductor?

I can't comment yet but from what you said it seems that you need to install additional libraries like the one you mentioned: "curl".

from the R terminal you should type:

install.packages("curl")

then R will install the package. Not sure if there are some other dependencies. If yes R will return you some error message in which there are some package names you need to install.

EDIT:

of course, after installing all the required packages you should be able to install biocLite too.


Need help with Pairwise Alignment module in iterating over the alignment from Bio import SeqIO
from Bio import Align

ref_seq_1 = SeqIO.read('C:/Users/King_Arthur/Desktop/ref_seq/segment 1/ref_seq_8.fasta','fasta')

seq1 = SeqIO.read('C:/Users/King_Arthur/Desktop/file/segment 1/Myfile_1 (1).fasta','fasta')

aligner = Align.PairwiseAligner()
aligner.mode = 'global'
aligner.match_score = 1
aligner.mismatch_score = -2
aligner.gap_score = -2

alignments = aligner.score(ref_seq_1.seq , seq1.seq)

print(alignments)
for alignment in sorted(alignments):
    print(alignment)

So this is my code and as you can see in the last section i am trying to iterate over my alignment but I am getting this error
TypeError: 'float' object is not iterable

I have tried various things like using str() but it gives some strange values and I also tried to read the source code by using the inspect module but I can't figure out the problem.
Any help would be really appreciated.
My final objective is to find out how many matches, mismatches and gaps are present in the final alignment using biopython.
if there is any other better way to do it in python please feel free to suggest.

aligner.score(ref_seq_1.seq , seq1.seq) doesn't produce an alignment. It just gives you the score of the alignment as float value.
Instead you need to run:
alignments = aligner.align(ref_seq_1.seq , seq1.seq)
Now alignments is an Bio.Align.PairwiseAlignments you can iterate over alignments or get the score with alignments.score which gives you the same as aligner.score(ref_seq_1.seq , seq1.seq)`
If you need more information, you could check the documentation.


Why Do I Get Different Results Each Time From AutoDock Vina? I just tried to use Vina for calculating binding affinity.
I used the same protein, the same ligand, the same options, but there are different results.
Why do I get different results?
Is it the limit Point of docking tools?
Then why we use it?

Vina uses a non-deterministic docking algorithm. See the manual.

The docking algorithm is non-deterministic. Even though with this receptor-ligand pair, the minimum of the scoring function corresponds to the correct conformation, the docking algorithm sometimes fails to find it. Try several times and see for yourself. Note that the probability of failing to find the mininum may be different with a different system.

The manual states that you can set the random seed in order to obtain reproducibility.

I changed something, and now the docking results are different. Why?
Firstly, had you not changed anything, some results could have been different anyway, due to the non-deterministic nature of the search algorithm. Exact reproducibility can be assured by supplying the same random seed to both calculations, but only if all other inputs and parameters are the same as well. Even minor changes to the input can have an effect similar to a new random seed. What does make sense discussing are the statistical properties of the calculations: e.g. "with the new protonation state, Vina is much less likely to find the correct docked conformation".


Outputting crunchable list of HBonds from Pymol X-posted from the Main Stack Overflow because folks here are possibly more familiar with PyMol(?)

I have been trying to assess the strength of an interface I want to mutate using Pymol. What I am interested in is a list of contacts withing a specified interface (which I have calculated using the InterfaceResidues.py script from the PyMol wiki script library). My aim is to produce a bunch of mutants using the mutagenesis wizard and assess how the HBond network changes.

Actions>Find>Polar Contacts does the trick in visualising them, but doesn't output a file with numbers which I can then crunch to assess which mutations weaken my interface most.
I have been trying the scripts list_contacts.py and list_hb.py from Prof. Robert Campbell of Queen's University at Kingston, but with limited success (nothing gets printed or written to the file).
The scripts can be found here http://pldserver1.biochem.queensu.ca/~rlc/work/pymol/

I am going to contact Prof. Campbell directly re: scripts but, in the meantime, and in case he's too busy to deal with a Pymol noob, does anyone have any ideas of what else I could try?

I am using Pymol 2.3.1 and Python 3.7.1

Thanks a lot to all who commented/replied, the author of the scripts provided a response. It was down to a syntax change in PyMol of which I wasn't aware that caused my selection to come up empty.


Different outputs of computeMatrix depending on binsize I'd like to output a figure in +/-1kb around a center using computeMatrix as following command:

computeMatrix reference-point -S data.bw -R centers.bed -a 1000 -b 1000 -bs=10 -out newdata10.gz --referencePoint center

The bed file is like:

chr1 100 101

If I change -bs=1, two outputs of the same input newdata10.gz are quite different. Which's correct with my commands? Maybe I don't understand the meaning of binsize option.

Please don't send things like this to the deepTools mailing list and post them here.

As I mentioned on the list, what you're likely being surprised by is that computeMatrix is not summing values, but taking the average in each bin. You can change this with the --averageTypeBins option (the default is mean, you presumably want sum).

Update: For those wondering the cause of the confusion was due to not using --missingDataAsZero. In other words, including or not including NAs will lead to different results, since different numbers of values will be included.


Plotting coverage of annotation over collection of region I'm trying to plot "meta" coverage of annotation: i.e. features (eg. gene class) over certain regions. It is similar to read coverage plots over gene body, except my input is two bed files (both in BED6 format) - (A) one containing the regions for which the meta plot is needed, and another being the list of genes or any feature whose coverage over the regions in (A) is needed.

Is there any package or tool which can create such plots (my domain is limited to python but I can try to work with R) ?

Something akin to this but for whole gene body[1]:

[1]- Kelley, D, Rinn, J (2012). Transposable elements reveal a stem cell-specific class of long noncoding RNAs. Genome Biol., 13, 11:R107.

Edit: The BED file containing list of regions looks like this:

chr3    39218734    40053659    region1 0.92426187419769    +
chr4    140163762   140453127   region2 0.896103896103896   -
chr7    40549151    41205036    region3 0.986072423398329   +
chr8    81291743    81963246    region4 0.94184168012924    -
chr9    12284032    12539789    region5 0.95539033457249    -

And the bed file containing features to be plotted looks like this:

chr3    39218100    40053200    LINE    1   +
chr4    140163962   140453027   LINE 1  -
chr7    40549002    41204999    SINE 1  +
chr8    81291143    81963846    LTR 1   -
chr9    12284332    12539720    LTR 1   -

If one assumes that repeat elements of a given type (e.g., LINEs) don't overlap each other, then the following will work:

Split your BED file by repeat element, such that you have a LINE.bed, SINE.bed, etc.
Convert those to bedGraph (e.g., awk 'BEGIN{OFS="\t"}{print $1,$2,$3,"1.0"}' LINE.bed > LINE.bedGraph).
Use UCSC tools to convert those bedGraph files to BigWig.
Install deepTools and run computeMatrix reference-point -b 2500 -a 20 -S LINE.bigWig SINE.bigWig LTR.bigWig -R Regions_Of_Interest.BED -o foo.mat.gz
Make a profile plot with plotProfile (plotProfile -m foo.mat.gz -o foo.png --perGroup)

If you do have overlapping regions you'll need to first make a disjoint set of intervals (you can use bedops for this), find the coverage of them (e.g., bedtools coverage) and then continue on with that.


how to estimate phylogeny without DNA sequences? I am new to bioinformatics, I am reading ‘Analysis of Phylogenetics Second Edition and Evolution with R’ from Emmanuel Paradis.
I can create phylogeny from DNA sequences, by first calculating the distances and then using (ie nj() function).
My question is: if instead of DNA I have trait data (ie height and weight of species), how can I estimate a phylogeny? (is the same process valid? somehow calculate distances and pass that to nj()?)

Sequence data is basically always better and at this point in history it is de rigueur for any serious study, but originally phylogenies were all based on (ideally discrete) phenotypic characters, so there's no reason why it can't work in principle. As @Pallie mentions in a comment, the clustering methods in e.g. NJ or UPGMA will work on any set of distance measures, and they were in fact originally developed on phenotypic characters rather than sequence data.
I'd recommend taking a look at chapter 24 of Joe Felsenstein's book on phylogenies, which is digitized here (yes the link looks sketchy, sorry). There is a section on p. 425 "Inferring phylogenies and correlations", which discusses how you can in principle use quantitative characters to infer phylogenies. The theory and math for why this is supposed to work is earlier in the chapter.
In many cases people prefer to code the characters as discrete because it makes it easier to work with. But there are also methods that work with the raw quantitative characters.
In R, you should be able to do something like this:
# some simulated data, species are rows and characters are columns
> data
      char1 char2 char3
spec1     0     3     6
spec2     1     4     7
spec3     2     5     8

> library(ape)
> tr = nj(dist(data))
> plot(tr)


Spearman correlation for large dataset I have two datasets (DataA and DataB) and I want to find the Spearman correlation between genes and also pull out the gene names (stored in first column of dataset) in R. I am using fread from read.table to read the file and cor.test to find Rho and p-value.
The code works fine for a smaller data set of 10 obs, 10 variables. But, my actual data set is larger (DataA - 795 obj, 542 variables and DataB - 925 obs, 542 variables).
When I try to run the code on my actual data, it is consuming a large amount of time (it has been 23 hrs since my code is running). Is there a way to optimize this code? Or is there a bug in the code?

library(data.table)
A.data <- fread("DataA.csv", header = TRUE, data.table = FALSE, stringsAsFactors = FALSE)
B.data <- fread("DataB.csv", header = TRUE, data.table = FALSE, stringsAsFactors = FALSE)

correlation.result <- data.frame()

## Computing Correlation

for (i in 1:nrow(A.data)) {
  for (j in 1:nrow(B.data)) {
       correln <- cor.test(as.numeric(A.data[i,2:ncol(A.data)]), as.numeric(B.data[j,2:ncol(B.data)]), method = "spearman", exact = FALSE)
       rho <- correln$estimate
       p.val <- correln$p.value
       Gene <- paste(A.data[i,1], B.data[j,1], sep="_", collapse=" ") # To pullout gene name pairwise
       comb.data <- data.frame(rho, p.val, Gene)
       correlation.result <- rbind(correlation.result,comb.data)
 }
}

Your code looks fine, but rbind to progressively build and resize a data frame can be gluttonous in terms of memory and time, especially given the A-row x B-row times you are doing this operation, which is fine for very small datasets, but not often realistic for real-world data.

Instead, perhaps cat results to a text file, where results are written in tab-delimited columns. You can then read in that text file later on.

In addition to serializing or saving your results to text, which make them easier to work with in non-R contexts, you won't have to deal with R's memory allocation problems when running rbind repeatedly. You just write out the answer to a file and read it in later on, in one pass.

Another approach is to preallocate memory to a presized data frame, and then write to that allocated space.

You know the size of your inputs ahead of time, which makes this approach tenable. In this case, it looks like you are storing two floats and a string (rho, p.val, and Gene) and that you are going to store A-row x B-row number of them.

Here's one way you might go about this:

# create a presized data frame
Arows <- nrow(A.data)
Brows <- nrow(B.data)
correlation.result <- data.frame(t(matrix(c(0.0, 0.0, "foo"), nrow=3, ncol=Arows * Brows)), stringsAsFactors=FALSE)
colnames(correlation.result) <- c("rho", "p.val", "Gene")

# assign values to that presized data frame
Acols <- ncol(A.data)
Bcols <- ncol(B.data)
for (i in 1:Arows) {
  # get A data
  Ai <- as.numeric(A.data[i, 2:Acols])
  Aname <- A.data[i, 1]
  for (j in 1:Brows) {
    # calculate row offset in presized data frame
    ij <- (i - 1) * Brows + j
    # get B data
    Bj <- as.numeric(B.data[j, 2:Bcols])
    Bname <- B.data[j, 1]
    # run test
    correln <- cor.test(Ai, Bj, method = "spearman", exact = FALSE)
    rho <- correln$estimate
    p.val <- correln$p.value
    Gene <- paste(Aname, Bname, sep="_", collapse=" ")
    # assign value to row of result data frame
    correlation.result[ij,] <- c(rho, p.val, Gene)
  }
}

This avoids R needing to resize and recopy the data frame contents as they grow, which I suspect is the major reason it is taking so long to run.

Chapter or "Circle" 2 of Patrick Burns' book The R Inferno discusses the problems with growing objects with rbind and ways to handle it — you might find this a useful read. This book is a PDF download via your search engine of choice.


You are looking at the supplementary data of a paper. That seems to have given you a list of features, and some information about those features. Specifically, you seem to have a list of two types of element:

"peg". Based on the information there, I assume "peg" stands for "protein encoding gene". Note that all peg lines have a peptide sequence (not seen in your question but present in the original link) and their "function" column is always some sort of protein function (e.g. "hypothetical protein" or "Transposase").

"repeat". These are presumably repetitive elements.

For each line, you are given an identifier for the sequence in which this element can be found. For example, your first line is on gi|1045318032|gb|MAZB01000333.1|. That is a GI (gi|1045318032) and a GenBank identifier (gb|MAZB01000333.1). If you go to the NCBI website, and search for 1045318032, you will find that it is the Candidatus Erwinia dacicola isolate Erw_SC contig_333, whole genome shotgun sequence.
So, your first line is telling you that the authors have identified a protein encoding gene on sequence gi|1045318032 and the gene starts at position 511 of that sequence and ends at position 2. This might sound odd, but that's because the gene is on the negative (-) strand and positions are calculated with respect to the positive strand. So, if you were to count on the negative, the gene would start at 2 and end at 511.
The same goes for all the other lines you have there. Just look up the GI number and you can know what sequences they are referring to.
Finally, the second column, the "feature_id" looks like something specific to this figshare service that is hosting the data. I dont' know what it is, but I would guess it's an identifier for this particular feature. I don't think it has any relevance outside http://figshare.com.


CCP4 file to a Python 3 numpy array or similar workaround I would like to merge together several ccp4 formatted density maps (and do a few minor things).

So ideally I would like to open the ccp4/mrc files as numpy arrays in Python 3 and save the array as a ccp4 file. CCTBX can open ccp4 maps as numpy arrays, but I really would like to avoid using Python 2.7 for various reasons.

This is to be used for constraining energy minimisation in PyRosetta using the pyrosetta.rosetta.protocols.electron_density.SetupForDensityScoringMover mover/sampler. But its pyrosetta.rosetta.core.scoring.electron_density.getDensityMap data loader only likes ccp4/mrc files (density maps) and not mtz files (structure factor maps).

From what I can tell there is no way to fudge a ccp4 map from a mtz map, correct? The fact that only ePDB provides ccp4 maps makes me want to use mtz files.

CCP4's Blend can merge datasets but uses only mtz files...

You can use Gemmi.

import gemmi
ccp4_map = gemmi.read_ccp4_map('my.map')
ccp4_map.setup()  # optional
np_array = numpy.array(ccp4_map.grid, copy=False)
# ...
ccp4_map.write_ccp4_map('out.ccp4')

The setup() call expands data to the whole unit cell.
The numpy array above accesses the data through Python's Buffer Protocol (it does not copy the data). So changing values in numpy array changes gemmi.Ccp4Map.

You can also use gemmi to read mtz file and transform map coefficients to a real-space map, see https://gemmi.readthedocs.io/en/latest/hkl.html.


How does the BWA-MEM algorithm assign its mapping qualities? Is there any resource (paper, blogpost, Github gist, etc.) describing the BWA-MEM algorithm for assigning mapping qualities? I vaguely remember that I have somewhere seen a formula for SE reads, which looked like

$C * (s_1 - s_2) / s_1,$

where $s_1$ and $s_2$ denoted the alignment scores of two best alignments and C was some constant.

I believe that a reimplementation of this algorithm in some scripting language could be very useful for the bioinfo community. For instance, I sometimes test various mapping methods and some of them tend to find good alignments, but fail in assigning appropriate qualities. Therefore, I would like to re-assign all the mapping qualities in a SAM file with the BWA-MEM algorithm.

Btw. This algorithm must already have been implemented outside BWA, see the BWA-MEM paper:

  GEM does not compute mapping quality. Its
  mapping quality is estimated with a BWA-like algorithm with suboptimal
  alignments available.

Unfortunately, the BWA-MEM paper repo contains only the resulting .eval files.

Update: The question is not about the algorithm for computing alignment scores. Mapping qualities and alignment scores are two different things:

Alignment score quantifies the similarity between two sequences (e.g., a read and a reference sequence)
Mapping quality (MAQ) quantifies the probability that a read is aligned to a wrong position.

Even alignments with high scores can have a very low mapping quality.

Yes, there bwa-mem was published as a preprint

  BWA-MEM’s seed extension differs from the standard seed extension in two aspects. Firstly, suppose at a certain extension step we come to reference
  position x with the best extension score achieved at query position y.

  ...

  Secondly, while extending a seed, BWA-MEM tries to keep track of the
  best extension score reaching the end of the query sequence

And there is a description of the scoring algorithm directly in the source code of bwa-mem (lines 22 - 44), but maybe the only solution is really to go though the source code.


Interpretation of hierarchical clustering with bootstrapping I have data that includes 'cases' and 'controls' and have carried out hierarchical clustering.

I used the pvclust package in R to bootstrap the results and significant branches are highlighted with red rectangles (based on au>0.95):

What is clear is that no clustering occurs that separates 'cases' and 'controls', and this is in fact what we expected and want to show. We want to show that the measured variable does not distinguish between cases and controls.

List itemApart from saying visually no clear clusters emerge that distinguishes between cases and controls are there any objective measures that can be used to say no significant clustering occurs between two groups?
One observation I have here is that the AU values and the BP values are very different, even though both p-values should be interpreted in a similar fashion, am I missing something?
Perhaps pvclust (bootstrapping) is not the right option here, is there a better way of showing quantifiably that no significant clustering occurs between two groups? Perhaps some kind of supervised clustering (I am not sure what this even means in this context)?

Good question. You have performed the first step via unsupervised learning and AU bootstraps give great results, orthodox bootstraps give no significant (or very few) clusters. The contrast between these approaches is unusual because they should be somewhere close to each other. If you can resolve that, your analysis is fine. So at present IMO the unsupervised learning needs further investigation.
... deleted

... antibodies ... oh I see. Trees, like you've done, don't work here because of the cross-reactivity creates network phylogenies. This is the antonym of bifurcating trees. If someone has a robust bifurcating phylogeny with antibody data and gets good bootstraps .. utterly amazing. However, if you subject your data to a network analysis you'll find connections all over the place - in practical terms this really messes up bootstraps. Thats the likely reason for you results. There is no point doing a complex bifurcating tree analysis because the antibody cross-reactivity will mess it up. k-means has been implemented as unsupervised learning in these instances and would provide the starting point in using clustering to define a supervised learning problem.


You see negative values with your function because you're setting the average of each row to 0 and its standard deviation to 1.

In general, I would trust a standard normalization method (rma in this case) more than some random "truncate and then scale the rows" method. Your method isn't even doing any between-array normalization, which is the benefit of rma.


Can we use non base-called fast5 files in poretools? I run the MinION MinKNOW without the live base-calling option. We know there is Metrichor and Albacore to perform base-calling after this process. However, I have not done any base-calling yet.

My question is: Is it possible to use directly the fast5 files with poretools to extract a fastq file without doing any previous base-calling process? I tried it and I get an empty fastq file.

My reads directory only contains the fast5 folder and I run:

poretools fastq fast5/ > output.fastq

Any ideas why do I get an empty file? What is the difference between base-called fast5 and non-base-called fast5 files?

Now I am trying to do base-calling with Albacore to see if I get a fastq file.

No, poretools does not do basecalling. The poretools fastq command can be used to extract the FASTQ information from the basecalled FAST5 file (via MinKNOW live-basecalling or albacore). Alternatively, both of these basecallers can export a FASTQ file directly if desired.

The difference between a basecalled FAST5 and a nonbasecalled FAST5 is the presence of the FASTQ within the file. This can be directly viewed using an HDF5 viewer, like h5dump from the commandline, or HDFView from the GUI.


Moving UMI from tag into read name I have a bam file with Unique Molecular Identifiers (UMIs) for each read present on the RX tag ( such as RX:Z:TGAGAAGGG), as expected by picard and fgbio tools.
However, several other tools (such as UMI-tools) expect the UMIs to be present at the end of the read name, such as (NB501943:73:HK333BGXC:3:12402:20313:16437_TGAGAAGGG).
Does anyone know of a tool to do this?
I have an awk solution which I'll self-answer with, but I don't know whether that is robust.

My current solution is (awk is not a tool I'm very familiar with, so this may not be good awk code, there is presumably a way to get it to print the header lines unmodified.):
samtools view -H temp.sam > temp2.sam
samtools view temp.sam |
awk '{OFS = "\t"} {for(i=1;i<=NF;i++)  if ($i ~ /^RX:Z:/) {$1=$1"_"$i; gsub("RX:Z:","",$1); print}}' >> temp2.sam

This seems to work when I test it, but I don't know if I'm missing anything.


Coronavirus phylogeny and evolution This new article presents the phylogeny of the coronaviruses and the placement of the new coronavirus in it. This made me wonder of what we know about the coronavirus evolution, particularly where it comes to evading human immune defenses: it is not mutating as fast as HIV - indeed, the genomes responsible for outbreaks in various countries are almost identical. In the same time it is unlikely to undergo gene reshuffling typical of influenza, since coronaviruses have single RNA genome.

FYI coronaviruses have proof reading activity, unlike HIV, its a huge genome, unlike HIV. Thus the point mutation rate is unlikely to be as rapid as HIV.

The phylogenies in this paper are not really about the evolution or epidemiology of the virus. The trees used for viral taxonomy and form part of the International Committee for the Taxonomy of Viruses (ICTV) assessment in officially naming the virus. The ICTV is the authority in naming a virus. The ICTV essentially called it SARS2 as opposed to the WHO who used the name COVID-19. Thats the point of this paper and will be massively cited.

The broader evolutionary epidemiology of "SARS2" is for further studies. The Lancet paper on 22nd February gave some indication of this (and I suspect you know about that paper).


Reduce number of transcripts in a highly variable de novo transcriptome assembly I have a de novo assembly using both multiple SRA and locally sequenced transcriptomes. I started with 270M PE reads from 9 tissues. Here are the assembly stats generated with TrinityStats.pl:

################################
## Counts of transcripts, etc.
################################
Total trinity 'genes':  543668
Total trinity transcripts:      1111089
Percent GC: 33.77

########################################
Stats based on ALL transcript contigs:
########################################

        Contig N10: 2117
        Contig N20: 1324
        Contig N30: 945
        Contig N40: 714
        Contig N50: 561

        Median contig length: 344
        Average contig: 496.17
        Total assembled bases: 551293315

#####################################################
## Stats based on ONLY LONGEST ISOFORM per 'GENE':
#####################################################

        Contig N10: 2172
        Contig N20: 1333
        Contig N30: 925
        Contig N40: 681
        Contig N50: 523

        Median contig length: 323
        Average contig: 473.63
        Total assembled bases: 257499729

As you can see, the total Trinity 'genes' number is very high. I suspect that this is due to a high biological variability of the transcriptomes (coming from organisms sampled in Italy, China, Spain...).

How can I assess whether the problem is really the biological variability?
How can I reduce this number? I'd rather not downsample the reads. It looks like CD-Hit might do what I'm looking for...

As you've suggested, CD-HIT works for reducing transcript numbers. We used a mixture of expression-based filtering and CD-HIT for reducing transcript counts for our genome-guided transcriptome assembly. This reduced numbers by a lot, without much change in BUSCO scores:

Map RNA-seq reads to Trinity-generated transcripts using Salmon
Use the expression of BUSCO genes to set a credible signal cutoff
Subset transcripts based on this threshold [for us it was 50 counts]
Identify the longest Met to Stop open reading frame for each transcript (may not be a good idea in all cases)
Run cdhit to collapse similar transcripts

Here's our CD-HIT command line:

cdhit -T 10 -c 0.98 -i longest_MetStopORF_HC50_TBNOCFED.fasta -o cdhit_0.98_LMOHC50_TBNOCFED.prot.fasta

More details here:

https://bmcbiol.biomedcentral.com/articles/10.1186/s12915-017-0473-4#Sec9


Will using smaller kmers help get larger contigs? If not, then what? I've been using MEGAHIT to assemble metagenomes, with particular focus on specific genomic areas.

Sometimes all I get is gene fragments or pathway fragments (eg. if I know that genes A, B, C, D and E should be together, I only get A, B, C in one contig, and maybe D and E in another). That is understandable when the sequencing depth is low, and I'm assume that it's the best I can get.

However. In two of my metagenomes, the genes of interest seem to be abundant (depth of 100-300+ listed for MEGAHIT contigs). And I still have the same problem. I assume that these 'breaks' between contigs are due to natural sequence variation which breaks the assembly process. (while sequencing errors are also possible, I've used bbduk to clean my data prior to assembly)

I'd like to have larger contigs because I'm curious about the gene order. Things I've tried (did not work): visualizing assembly graphs with Bandage around Blast hits, starting with a smaller k-min (21 as opposed to 27).

So I guess I'm wondering if using even smaller k-mers might make a difference. Or if someone has another suggestion, I'm happy to listen! Thank you.

Megahit setting used:

megahit -r input_files.fastq --num-cpu-threads 32 --min-contig-len 300 --presets meta-large -o output
megahit -r input_files.fastq --num-cpu-threads 32 --min-contig-len 1500 --presets meta-sensitive -o output
(meta-large starts from minimum kmer size 27; meta-sensitive from 21)

You can try s-aligner (free for 15 days) which usually gets quite larger contigs and quite larger NG50 for metagenomic data containing viruses (also phages). Indeed metaSPAdes also gets larger contigs, despite shorter than s-aligner.

Ref. s-aligner: a greedy algorithm for non-greedy de novo genome assembly


GROMACS Coordinates and velocity and force trajectories do not match accuracy I am using GROMACS to perform MD simulations of proteins. So, I tried to record the simulation results nstxout, nstvout, and nstfout with the same step size (nstxout 1, nstvout 1, nstfout 1). Then, I calculated the coordinate (step n) of an arbitrary step of the trajectory (step n) from the velocity of step n - 1, step n, and calculated the velocity of n - 1/2, and used that velocity to calculate the coordinate (n-1) of step n using the velocity of n - 1/2. LeapFrog was used in the simulation.
How can I get the coordinates, velocity, and force to be consistent?

From GROMACS documentation about reproducibility:

The following factors affect the reproducibility of a simulation, and
thus its output:

Precision (mixed / double) with double giving “better” reproducibility.
Number of cores, due to different order in which forces are accumulated. For instance (a+b)+c is not necessarily binary identical
to a+(b+c) in floating-point arithmetic.
Type of processors. Even within the same processor family there can be slight differences.
Optimization level when compiling.
Optimizations at run time: e.g. the FFTW library that is typically used for fast Fourier transforms determines at startup which version
of their algorithms is fastest, and uses that for the remainder of the
calculations. Since the speed estimate is not deterministic, the
results may vary from run to run.
Dynamic linking to different versions of shared libraries (e.g. for FFTs)
Dynamic load balancing, since particles are redistributed to processors based on elapsed wallclock time, which will lead to (a+b)+c
!= a+(b+c) issues as above
Number of PME-only ranks (for parallel PME simulations)
MPI reductions typically do not guarantee the order of the operations, and so the absence of associativity for floating-point
arithmetic means the result of a reduction depends on the order
actually chosen
On GPUs, the reduction of e.g. non-bonded forces has a non-deterministic summation order, so any fast implementation is
non-reprodudible by design.

I would add that in your case of hard coding everything by hand you may have to look a bit into your barostat and thermostat (although their influence should be found on speed of the atoms). Maybe a sample code with coordinates and speed of 4 atoms would help to give a more precise answer.
You can also try having a look at gromacs source code. I believe OpenMM provides a bit more control and clarity on each step, if for some unfathomable reason you are not doing this only for curiosity sake but because you actually have to make a code that reproduces these results.


Linear models of complex diseases A popular framework to analyze differences between groups, either experiments or diseases, in transcriptomics is using linear models (limma is a popular choice).

For instance we have a disease D with three stages as defined by clinicians, A, B and C. 10 samples each stage and the healthy H to compare with is RNA-sequenced. A typical linear model would be to observe the three stages~A+B+C independently. The data of each stage is not from the same person. (but for the question assume it isn't)

My understanding is that such a model would not take into account that stage C appears only on 30% of patients in stage B. And that a healthy patient upon external factors can jump to stage B.

If we want to find the role of a gene in the disease we should include somehow this information in the model. Which makes me think about mixing linear models and hidden Markov chains.

How can such a disease be described in terms of linear models with such data and information?

There are two potential sources of bias in this design.

We cannot distinguish correlation from causation.

Imagine two cases. In the first, the disease progression is inducing immune response. Later stages will be associated with the higher gene expression levels. In the second scenario, the disease is caused by overexpression of a gene. Later stages will be also associated with the higher expression.

This is typical for observational studies. But I just want to mention that a special care should be taken during interpretation of the results.

If we are not following our individuals, we cannot distinguish correlation and avoidance.

Let's say the disease is lethal in certain cases; the survival is negatively correlated with the disease stage. Now imagine there is a gene, which causes severe symptoms when highly expressed. On the later stages you will only observe those patients in which the gene was not highly expressed. From that you would conclude that gene expression is decreasing with the disease progression. In reality this gene is very important and causal, you just do not have patients which are alive to see this.

This is similar to Wald's studies of aircrafts (Survival bias).

  Researchers from the Center for Naval Analyses had conducted a study
  of the damage done to aircraft that had returned from missions, and
  had recommended that armor be added to the areas that showed the most
  damage.
  Wald proposed that the Navy instead reinforce the areas where the
  returning aircraft were unscathed, since those were the areas that, if
  hit, would cause the plane to be lost.

I think the second point is crucial, and can and will lead to false conclusions. I do not see immediately how to avoid this bias if you study different individuals.

Therefore I suggest that the same individuals are followed for a long time.

There are different approaches you can use later. For example you can have two-stage procedure:

Identify genes which are differentially expressed between healthy (H) and sick (A, B or C).
Build a linear model of disease stage stage ~ gene1 + gene2 + ..., using genes identified at step 1.
Similarly build a linear model of survival as a function gene expression.
It is possible you can use logistic regression to infer probability of transition from stage B to stage C.

This is not the only possible approach of course. As you suggested, Markov model is also applicable.

While following the same individuals it is possible to use continuous-time Markov model for such a study. In this case discrete or continuous measurements are recorded at certain time points, and the model parameters are inferred using maximum likelihood.

I'm not an expert in this field, but I think this paper describing the msm package for R will be useful. It also supports hidden Markov models.


How to fragment genomes into non-overlapping sequences of differing sizes? I have bacterial chromosomes/plasmids that I have downloaded off refseq in a multi-fasta that I would like to fragment into non-overlapping fragments of between 1 kb and 15 kb.
I have tried googling for related tools but to no avail. Does anyone know how I would achieve this? Something very similar is done in this paper:
https://link.springer.com/article/10.1007/s40484-019-0187-4
Thanks in advance!

You could use numpy.random.randint to sample integers from a uniform distribution, from 1000 to 15000 nt, extracting subsequences of that length from a linearized multi-FASTA file.
Here's some code:
#!/usr/bin/env python

import sys
import numpy as np

low = int(sys.argv[1])
high = int(sys.argv[2])

def process_record(header, sequence, low, high):
    next = 0
    fragment = 1
    while next < len(sequence):
        new_next = np.random.randint(low, high) + next
        if new_next > len(sequence) - high:
            new_next = len(sequence)
        sys.stdout.write('{}:{}\n{}\n'.format(header, fragment, sequence[next:new_next]))
        next = new_next
        fragment += 1

header = ""
sequence = ""
for line in sys.stdin:
    if line.startswith('>'):
        if len(sequence) > 0:
            process_record(header, sequence, low, high)
            sequence = ""
        header = line.rstrip()
    else:
        sequence += line.rstrip()

if len(sequence) > 0:
    process_record(header, sequence, low, high)

To use it:
$ ./13879.py 1000 15000 < in.fa > out.fa

Fragments will be disjoint (non-overlapping, contiguous) subsequences.
I add a check to determine if the last fragment of a FASTA record would be shorter than low (1kb) in length, and, if so, append that fragment to the final subsequence. It is possible for that last fragment to be up to 16 kb - 1 nt in length, given your parameters:
if new_next > len(sequence) - high:
    new_next = len(sequence)

You could change this logic to discard the last fragment (i.e. break), or simply allow fragments less than 1kb in size, depending on what you need.
EDIT
Here's a slight modification that should allow to you preserve fragment sizes within your bounds.
A cartoon may help illustrate the problem. Let's set up a FASTA record, which has been broken up into n fragments:
   1       2           n-2      n-1      n
|-----|---------|...|-------|---------|-----|

We may get unlucky and sample a random integer for fragment n-1 (frag(n-1)) which makes the length of frag(n) less than low (1 kb in your example).
In my first suggested solution, frag(n-1) and frag(n) are appended to one another. This is why you get some fragments longer than high (15 kb).
Another option is to discard frag(n), which is of length less than low (1 kb). Here's some code:
#!/usr/bin/env python

import sys
import numpy as np
import timeit

low = int(sys.argv[1])
high = int(sys.argv[2])
in_fn = sys.argv[3]

"""
Globals
"""
fragment_lengths = []

"""
Process a record
"""
def process_record(header, sequence, low, high):
    next = 0
    last = 0
    fragment = 1
    fragment_sequence = ""
    while last < len(sequence) - high:
        next = np.random.randint(low, high + 1) + last
        fragment_sequence += sequence[last:next]
        # Comment the sys.stdout.write and uncomment the fragment_lengths line, if measuring fragment length distributions
        sys.stdout.write('{}:{}-{}:{}\n{}\n'.format(header, last, next, fragment, sequence[last:next]))
        #fragment_lengths.append(next - last)
        last = next
        fragment += 1
    """
    Uncomment to compare input sequence with fragments
    """
    #sys.stderr.write('In > {}\nOut > {}\n'.format(sequence, fragment_sequence))

"""
Parse records
"""
def parse_records(in_fn, low, high):
    header = ""
    sequence = ""
    with open(in_fn, 'r') as ifh:
        for line in ifh:
            if line.startswith('>'):
                if len(sequence) > 0:
                    process_record(header, sequence, low, high)
                    sequence = ""
                header = line.rstrip()
            else:
                sequence += line.rstrip()
        if len(sequence) > 0:
            process_record(header, sequence, low, high)

def count_elements(seq):
    hist = {}
    for i in seq:
        hist[i] = hist.get(i, 0) + 1
    return hist

def wrapper(func, *args, **kwargs):
    def wrapped():
        return func(*args, **kwargs)
    return wrapped

def main(in_fn, low, high):
    """
    Comment below if testing the distribution of fragments
    """
    parse_records(in_fn, low, high)
    """
    Uncomment to test the distribution of fragments from
    running fns on same input 100k times
    """
    #wrapped = wrapper(parse_records, in_fn, low, high)
    #timeit.timeit(wrapped, number=100000)
    #print(count_elements(fragment_lengths))

if __name__ == "__main__":
    main(in_fn, low, high)

If frag(n) is needed, and you need to preserve the uniform distribution (or whatever fragment length distribution you want to model), then rejection sampling might be an option.
With rejection sampling, you keep fragmenting your input sequence, until one of the distributions of lengths gets you a set of fragments that entirely cover the input sequence.
If you need that, follow up and I'll take another look.


Filter VCF for "CAF=[0.9,0.1]" Related to this other question
but in this case the CAF is defined as "An ordered, comma delimited list of allele frequencies...".

I tried using vcffilterjdk and the command suggested by Pierre

java -jar dist/vcffilterjdk.jar -e 'variant.variant.getAttributeAsStringList("CAF","1.0").stream().skip(1L).map(A->Double.parseDouble(A)).anyMatch(A->A<0.01);' input.vcf

but predictably returned an error

java.lang.NumberFormatException: For input string: "0.1]"

I don't know javascript so how should I change the command to make it work?

EDIT:

HEADER

##INFO=<ID=CAF,Number=.,Type=String,Description="An ordered, comma delimited list of allele frequencies based on 1000Genomes, starting with the reference allele followed by alternate alleles as ordered in the ALT column. Where a 1000Genomes alternate allele is not in the dbSNPs alternate allele set, the allele is added to the ALT column.  The minor allele is the second largest value in the list, and was previuosly reported in VCF as the GMAF.  This is the GMAF reported on the RefSNP and EntrezSNP pages and VariationReporter">

LINE

chr1    51479   rs116400033     T       A       .       PASS    ASP;CAF=[0.8926,0.1074];COMMON=1;G5;G5A;KGPROD;KGPhase1;OTHERKG;RS=116400033;RSPOS=51479;SAO=0;SSR=0;VC=SNV;VLD;VP=0x050000000005170016000100;WGT=1;dbSNPBuildID=132

Printing a bracketed list of numerical values is not typical for VCF—I'd guess it's not strictly allowed. A list should be printed only with the commas to separate values. Perhaps the solution is as simple as removing the square brackets.

cat input.vcf | perl -ne 's/CAF=\[([^\]]+)\]/CAF=$1/; print' > fixed.vcf

Replace cat with gunzip -c if your VCF file is gzip-compressed.

The Perl regular expression is

looking for the string CAF=[...]
grabbing the characters in ...
replacing the entire match with CAF=...

Here's a detailed breakdown of the Perl regular expression.

s/.../.../; — substitution syntax; replace the contents of the left block with the contents of the right block
[^\]]+ — one or more non-] characters; this is the trickiest bit; see this tutorial for more info on the [] syntax for character classes and the special ^ symbol for character class negation
(...) — the parentheses capture the matched pattern into the $1 variable
CAF=$1 — replaces everything matched in the first block with CAF= followed by the captured characters


VCF merge containing CNV How do I merge VCFs files containing CNVs?

I use vcf-merge, a VCFtools function, and after bgzip and tabix, SAMtools, to index and tab separate variants, but I don't know if it is the right way.

vcf-merge (or better: bcftools merge) does a column-based merge of VCF files, basically taking on an additional genotype column to each variant. While this sort-of works, it doesn't take into account all the information about variation. There might be low coverage at a region in one VCF file that indicates a deletion, but the deletion is only observable when considering other samples with high coverage at that same region. It's also possible that INDELs could be defined differently in different samples, requiring an additional bcftools norm step for left-normalisation prior to the merge.

Assuming this column-based merge is what you want (and if not, then bcftools concat is your friend for row-based merging), then it's better to regenerate the merged VCF file as a multi-column file directly. One way of doing this is feeding all BAM files into samtools mpileup to generate genotype likelihoods, then feeding the result of that into bcftools call. Something like the following:

samtools mpileup -uf ref.fa in1.bam in2.bam in3.bam | \
  bcftools call -mv -O z -o variants_ref_1-3.vcf.gz


Associating SNP and GENE Assuming I have SNPs data using hg19, how can I know which SNP belongs into which Gene?

The data looks like:

chr10_103577643
chr10_124712463

and so on.

I want to add a column of Gene, which would tell to which Gene the SNP belongs.
The file is not a VCF file but a Txt file.


Is there a single cell RNAseq equivalent of GTEx or TCGA? Or do I need to find individual studies and obtain data the long way!

I think the most well-known scRNA dataset that is similar to GTEx in the Human Cell Atlas (https://www.humancellatlas.org). It is pretty much still ongoing but there are some preliminary data.


Comparing two files with overlapping regions and get their associated information I have two file that have overlapping regions. One file has chromosomes, position, CNV data and samples.Other file has chromosomes, positions and ensemble ids of genes.

Here is some lines of two files.

contig  start   end     type  bi-allelic  S1  S2  S3
chr1    12048   177033    DUP FALSE   -1  22  24
chr1    12048   89237     DUP FALSE   12  10  -1
chr1    2985712 3355300   DEL   TRUE  2   2   2

Second file is like this

chr1    901877  911245  gene_id "ENSG00000187583.6"
chr1    22013   41670 gene_id "ENSG00000160075.10"
chr1    2985732 3355185 gene_id "ENSG00000142611.12"

I want output like this.

chr1    22013   41670 gene_id "ENSG00000160075.10"  DUP FALSE   12  10  -1
chr1    2985732 3355185 gene_id "ENSG00000142611.12" DEL   TRUE  2   2   2

Is there any way to do this?


shiny code with seurat object I am trying to create a shiny board for my scRNA seq data. I have saved my seurat objects adn I use following code as below:
library(shiny)
library(shinythemes)
library(Seurat)
library(uwot)

set_1 <- readRDS("/temp2/data/mydataB10.rds")
set_2 <- readRDS("/temp2/data/mydataB5.rds")

ui <- fluidPage(

    # Application title
    titlePanel("scRNAseq data browser"),

    fluidRow(
        column(4,

               sidebarPanel(
                   radioButtons("dataset", label = h3("Dataset"),
                                choices = list("sample1" = "set_1", "sample2" = "set_2"),
                                selected = "set_1"),
                   helpText("Enter gene names in CAPS.  Gene names must be exact."),
                   textInput("gene2", label = "Gene Name", value = "GAPDH")

               )
        ),
        column(8,
               column(6,
                      mainPanel(
                          plotOutput("dimPlot1", width = "150%"),
                          plotOutput("genePlot2", width = "150%")
                      )
               )
        )
    )
)

server <- function(input, output, session) {

    datasetInput <- reactive({
        df <- (input$dataset)
            })

    output$dimPlot1 <- renderPlot({
        DimPlot(datasetInput(), reduction = "umap")
    })

    output$genePlot2 <- renderPlot({
        FeaturePlot(datasetInput(), features = (input$gene2), reduction = "umap")
    })

}

shinyApp(ui, server)

That error that I get on the section where image needs to be displayed is "Subscript out of bound".
Following dataset can be downloaded and used as an example curl -Lo pbmc_10k_v3.rds https://www.dropbox.com/s/3f3p5nxrn5b3y4y/pbmc_10k_v3.rds?dl=1
Can anyone help in fixing this shiny script so that it works.
Many thanks.

Not easy to reproduce your error but here is my guess. This is a string, you can check it:
input$dataset
class(input$dataset)

This is a simple solution. If you want to call out the dataset depending on input, use a list:
set_1 <- readRDS("/temp2/data/mydataB10.rds")
set_2 <- readRDS("/temp2/data/mydataB5.rds")
data_list = list(set_1=set_1,set_2=set_2)

Then in the server :
datasetInput <- reactive({
        df <- data_list[[input$dataset]]
            })


Does this FASTQ data contain single or paired end calls? I have this fastq data from GEO:

zcat SRR1658526.fastq.gz | head -n 20
@SRR1658526.1 HWI-ST398:296:C1MP4ACXX:1:1101:1093:2094 length=102
GATCTCTATTACTTTTTGAAGGATTNNNNNNNNNNAANTTTTGAATCANNNNNNNNNNNNNNNNNNNNNNNNNNNNNTNNNNNNNNNNNNNNNNNNNNNNNN
+SRR1658526.1 HWI-ST398:296:C1MP4ACXX:1:1101:1093:2094 length=102
@<@FFDFDHHFHHIIIIGHBGGGIG##########10#0:BGDDHHII######################################################
@SRR1658526.2 HWI-ST398:296:C1MP4ACXX:1:1101:1167:2107 length=102
ATAATATTGTAGATATAAATGTTATCTAATCTTATCTGATCAGCTTGCTNNATANNNNNNNNNNNNNNNNNNACNTATGNNNNNNNNNNNNNNNNNNNNNCC
+SRR1658526.2 HWI-ST398:296:C1MP4ACXX:1:1101:1167:2107 length=102
CCCFFFFFHHFHHJJIIIIJIJJHJHJJJJIJJIJJJJIJHIJJJJGII#####################################################

It is supposed to be paired-end sequences. Is the prefix @ and + are the R1 and R2? What's the convention here?

Entries in a fastq file occupy 4 lines each and R1 and R2 are typically in different files. Since that SRA project is 2x51 you seem to have run fastq-dump without the --split-3 option, so both R1 and R2 are merged together. Make your life easier and never use SRA, but instead just download the individual files from ENA.


Psipred installation I am trying to install psipred in a unix server.

Psipred is installed, but fails because it can not find the correct database files.  I have downloaded uniprod90.fasta. Psipred also loads the old legacy ncbi BLAST-2.2.26 and BLAST+.
runpsipred and runpsipredplus have been patched to call the correct version of BLAST.

I have run a few commands to index the uniprod90 database:

pfilt uniref90.fasta > uniref90filt
formatdb -t uniref90filt -i uniref90filt
makeblastdb -in uniref90.fasta -dbtype prot -input_type fasta -out uniref90.fasta

I then ran psipred and got this error:

$ ./psipred/example/example.fasta
Running PSI-BLAST with sequence  ./psipred/example/example.fasta ...
[blastpgp] WARNING: Unable to open uniref90.pin
[blastpgp] WARNING: Unable to open uniref90.pin
FATAL: Error whilst running blastpgp - script terminated!

makeblastdb created 26 pin files, but not uniref90.pin.

Try with:

makeblastdb -in uniref90.fasta -dbtype prot -input_type fasta -out uniref90

The -out is the name of the DB. So if you use uniref90.fasta, it creates uniref90.fasta.pinwhile the script seems to be expectinguniref90.pin`.


What is the distribution of indel sizes in a healthy human genome? of insertion:deletion ratios? My understanding is that indels are from 1bp to 10Kb, and a healthy genome has ~400K-500K Indels. Surely most of these are small.

What is the distribution of insertion sizes in a healthy human genome? What is the distribution of deletion sizes? What is an average ratio of insertions:deletions?

I have not found a definitive reference.

One of the 2015 papers from the 1000 genomes project has a nice figure (figure 1) showing the size distribution of medium to large sized insertions and deletions:

From another 2015 1000 genomes paper, one can see that the absolute number of smaller indels is much larger, though an exact size range isn't given (as far as I saw). If you really want to know that, just download the most recent 1000 genomes VCF file(s) and compute the exact median size and/or distribution that you want.


Receptor-ligand database How can I find ligands of receptors (or vice versa)?

For example, I want to find the receptors of IL-8 (CXCL8).

CellPhoneDB is a publicly available repository of curated receptors, ligands and their interactions.

The interactions can be searched or downloaded.


How to download gene expression data from NCBI gene database In the NCBI gene database, I can add the expression tracks (circled in picture blow) through 'Tracks' button, but How I can download the expression data directly, not just look the picture?

You can download raw data from SRA (NCBI). On this website under RNA-Seq alignments, you'll find the samples. Click on them, for example SAMN05231885, and then click on PRJNA325427 to see the link to the SRA data (click on 1). Here you can download the raw fastq files of the experiment.


DiffBind: Model must include count data for contrasts. I am trying to use DiffBind for my ATACseq experiment going through the vignette for DiffBind. I was able to import the data and make a heatmap with plot(data). But somehow calling data <- dba.contrast(data, categories=DBA_TISSUE) function is not working for me giving a warning message:

  Warning message:
  Model must include count data for contrasts.

I tried running beforehand: data <- dba.count(data) but it is giving an error:

  Error: Error processing one or more read files. Check warnings().
  In addition: Warning messages:
  1:
  2:
  3:
  4:
  5:
  6:
  7:

My sample_sheet.csvlooks like that:

I just need to compare Dura vs Microglia which is Tissue column. How could this error be fixed? Any suggestions would be greatly appreciated.

The problem was that I was using wrong files. I ran a pipeline from here and it generates several folders with many files. The .bam files that should be used are located in the aligned_mm10 folder and are with the extension: .pe.q10.sort.rmdup.bam.


Calculating Synonymous and non-synonymous library(seqinr)
library(ape)
library(phangorn)

alignment file data file
    sylvia.aln <- read.alignment("abhi_seq/all_seq_no_gal.fas", format="fasta")
    sylvia.aln

kaks(sylvia.aln, verbose=TRUE)

Warning message:
In kaks(sylvia.aln) : sequence lengths are not a multiple of 3

It ends with this warning and there are no result.

When i ran one of the seqinr example data what I see there are no gaps in the sequence files. If i understand the results is due to the gap I suppose.

Any help or suggestion would be really appreciated.

Woaa that is complicated data.

The reason it fails is easy, correcting it is harder. It fails because the open reading frame is lost during in silico translation, e.g. several gaps have 4 indels and the final nucleotides after the final indel are always 4. Thus it will invariably encounter a stop codon OR it will simply refuse to deal with an indel if it is incorporated into a triplet codon. There is no triplet codon model for an indel so the whole calculations stops.

Correcting data that complex ... the easiest way to do it is align you entire data set using amino acids, e.g. muscle. Have an identical nucleotide data set. MAKE SURE the names and positions are identical between your amino acid and nucleotide data, as a programmer this isn't hard to do however, everyone has their own scripts to do this. If you don't do this you next step will fall over.

Use a amino acid to nucleotide alignment translator. I use transalign (I think its EMBOSS). This will give you an inframe alignment. I noted many/all of you start codons are ATG which is a good sign. Check the nucleotide output and make sure there is no overhangs of the triplet codon at the terminal 3' end and check in Jalview or whatever your alignment editor of choice is that the translation is okay.

This will alignment work in the synonymous/non-synonymous analysis and all subsequent/complimentary analyses.

If you have an intron in protein your you'll need to remove it (hopefully not).

Finally you need to be careful of '3rd codon saturation' in your data because that heavily affects your trees and its certainly a risk given the heterogeneity of the indel structure, but first I would get an inframe alignment. There are checks for saturation (I'm sure even R can do this).


Error in quantile.default(x, p = p) in EgdeR calcNormFactors I am trying to run TMM normalization using rpy2 and when I run calcNormFactors() function:

dge_list = edgeRLib.DGEList(counts=rawcounts)
dge_list = edgeRLib.calcNormFactors(dge_list, method="TMM")

I am getting the error:

  RRuntimeError: Error in quantile.default(x, p = p) :
    missing values and NaN's not allowed if 'na.rm' is FALSE

So, I thought that maybe my rawcounts 2D numpy array contains nan values but it does not. I checked it with both pd.isnull() and np.isnan() functions. So, I suspect now that the 2D array itself can not be normalized for some reason and I do not know why. The error is most probably not related to rpy2 set up and R installation because it was working before for similar numpy arrays but failed this time. Any suggestions would be greatly appreciated.

Ok, it turned out that some rows and/or columns contained all values as zeros and that was the reason for the error. I wrote a function to fix that:

def fix_rawcounts(rawcounts, cell_ids, gene_symbols):
    '''The function checks whether there are
    any rows or columns with all zeros
    and removes them'''
    print('Checking whether rows or columns contain all zeros')

    # Checking columns
    count_arr = np.apply_along_axis(np.count_nonzero, axis=0, arr=rawcounts )
    rawcounts = rawcounts[:,count_arr != 0]

    # Checking rows
    count_arr = np.apply_along_axis(np.count_nonzero, axis=1, arr=rawcounts )
    rawcounts = rawcounts[count_arr != 0,:]

    print('finished checking zero rows and columns')
    return rawcounts


Recommended coverage depth for ChIP Seq (H3K9me3 ) According to ENCODE recommendations

For broad-peak histone experiments, each replicate should have 45 million usable fragments

According to my understanding, a usable fragment means a uniquely mapped read / pair of reads.
However

H3K9me3 is an exception as it is enriched in repetitive regions of the genome. ... Tissues and primary cells should have 45 million total mapped reads per replicate.

I fail to grasp Encode's language.

Does "total mapped reads" mean uniquely mapped reads?
If so, why is 45 million "an exception" ? It is the same as the other broad peaks.

Any ideas where else can I look for an explanation?

In the initial description for "usable fragments", ENCODE has the following (emphasis mine):

Usable fragment – A fragment is defined as the sequencing output corresponding to one location in the genome. If single-ended sequencing is performed, one read is considered a fragment. If paired-ended sequencing is performed, one pair of reads is considered a fragment. Fragments are considered usable if they pass the various filters in the ChIP-seq uniform processing pipelines. Used to evaluate ChIP-seq data.

Aligning reads to repetitive regions is difficult because they don't give unique alignments.
Reads from these regions will often have poor alignment scores and may not pass alignment QC checks.
Reads may also be counted as multi-mapped, and often methods call for uniquely mapped reads, which again, may not be true for a large number of reads from H3K9me3 ChIP-seq.
So for this particular histone mark, ENCODE loosens the definition for what you should count towards a minimum sequencing depth threshold.