text:
calculating nucleotide frequency per column
I have some sequences shown below
CAGGTAGCC
CCGGTCAGA
AGGGTTTGA
TTGGTGAGG
CAAGTATGA
ACTGTATGC
CTGGTAACC
TATGTACTG
GCTGTGAGA
CAGGTGGGC
TCAGTGAGA
GGGGTGAGT
TGGGTATGT
GAGGTGAGA
CAGGTGGAG

Each line has 9 nucleotides.
Consider it to be 9 columns.I want to calculate the nucleotide frequency of each nucleotide for each of the 9 columns.
For example 1st column will have these bases C,C,A,T,C,A etc
Out put should be something like this
A   0.25    0.34    0.56    0.43    0.00    0.90    0.45    0.34    0.31
C   0.45    0.40    0.90    0.00    0.40    0.90    0.30    0.25    0.2
G   0.00    0.00    0.34    1.00    0.30    0.30    0.35    0.90    0.1
T   0.24    0.56    0.00    0.00    1.00    0.34    0.45    0.35    0.36

Note, I just made up the numbers to show you the output file format

I think a package in Bioconductor can do this, I don't recall exactly one. Has you searched there? Or what have you tried to calculate this? The original file is a plain text or a fasta file ( I don't recall if it is only for fasta or for other file types)?

Original file is a text file only. I am thinking of using awk(associative arrays) to calculate it, still writing the script though

For anyone interested in doing this from any sort of alignment file, I've implemented a position frequency matrix function in AlignBuddy.
input:
 15 9
seq_1      CAGGTAGCC
seq_2      CCGGTCAGA
seq_3      AGGGTTTGA
seq_4      TTGGTGAGG
seq_5      CAAGTATGA
seq_6      ACTGTATGC
seq_7      CTGGTAACC
seq_8      TATGTACTG
seq_9      GCTGTGAGA
seq_10     CAGGTGGGC
seq_11     TCAGTGAGA
seq_12     GGGGTGAGT
seq_13     TGGGTATGT
seq_14     GAGGTGAGA
seq_15     CAGGTGGAG

Command:
:$ alignbuddy input.phy --pos_freq_mat

Output:
### Alignment 1 ###
A       0.133   0.400   0.133   0.000   0.000   0.400   0.467   0.067   0.400
C       0.400   0.267   0.000   0.000   0.000   0.067   0.067   0.133   0.267
G       0.200   0.200   0.667   1.000   0.000   0.467   0.200   0.733   0.200
T       0.267   0.133   0.200   0.000   1.000   0.067   0.267   0.067   0.133

awk
awk '{L=length($1);for(i=1;i<=L;i++) {B=substr($1,i,1);T[i][B]++;}} END{for(BI=0;BI<4;BI++) {B=(BI==0?"A":(BI==1?"C":(BI==2?"G":"T")));printf("%s",B); for(i in T) {tot=0.0;for(B2 in T[i]){tot+=T[i][B2];}printf("\t%0.2f",(T[i][B]/tot));} printf("\n");}}' input.txt

A   0.13    0.40    0.13    0.00    0.00    0.40    0.47    0.07    0.40
C   0.40    0.27    0.00    0.00    0.00    0.07    0.07    0.13    0.27
G   0.20    0.20    0.67    1.00    0.00    0.47    0.20    0.73    0.20
T   0.27    0.13    0.20    0.00    1.00    0.07    0.27    0.07    0.13

Or, expanded for clarity:
awk '{
        L=length($1);
        for(i=1;i<=L;i++) {
            B=substr($1,i,1);
            T[i][B]++;
        }
      }
      END{
            for(BI=0;BI<4;BI++) {
                B=(BI==0?"A":(BI==1?"C":(BI==2?"G":"T")));
                printf("%s",B);
                for(i in T) {
                    tot=0.0;
                    for(B2 in T[i]){
                        tot+=T[i][B2];
                    }
                printf("\t%0.2f",(T[i][B]/tot));
                }
            printf("\n");
            }
}' input.txt

Can we do same for taking dinucleotide at time (eg; AA, TT, CG etc); the output matrix would reduce to half?

Here is one example of how to do this with a bit of python. Alternatively one could create strings of each column and using letterFrequency() from the Biostrings package.
#Make a list of hashes
hl = []
for i in range(9):
    hl.append({'A': 0, 'C': 0, 'G': 0, 'T': 0})

f = open("foo.txt")  # CHANGE ME
nLines = 0
for line in f:
    for idx, c in enumerate(line.strip()):
        hl[idx][c] += 1
    nLines += 1
f.close()

nLines = float(nLines)
for char in ['A', 'C', 'G', 'T']:
    print("{}\t{}".format(char, "\t".join(["{:0.2f}".format(x[char]/nLines) for x in hl])))

The output of your example is then:
A   0.13    0.40    0.13    0.00    0.00    0.40    0.47    0.07    0.40
C   0.40    0.27    0.00    0.00    0.00    0.07    0.07    0.13    0.27
G   0.20    0.20    0.67    1.00    0.00    0.47    0.20    0.73    0.20
T   0.27    0.13    0.20    0.00    1.00    0.07    0.27    0.07    0.13

how to change enumerate if I want to take dinucleotides at time (eg; AA, TT, CG etc); the output matrix would reduce to half?

Here's a Perl approach:
#!/usr/bin/env perl
use strict;
my %counts;

## Read the input file line by line
while (my $line = <>) {
  print;
  ## remove trailing '\n' characters
  chomp $line;
  ## split the line into an array at every character
  my @columns=split(/./, $line);
  ## iterate over the array from the first to the last position
  for my $i (0..$#columns){
    ## The nucleotide at this position
    my $nt = $columns[$i];
    ## Save the count in this hash of arrays. The keys are the
    ## nucleotides, and the value at each position is increased
    ## when that nucleotide is found in that position.
    $counts{$nt}[$i]++;
  }
}
## Iterate over the counts hash
for my $nt (sort keys(%counts)){
  print "$nt\t";
  ## dereference the array stored in the hash
  my @countsForThisNt = @{$counts{$nt}};
  ## iterate over the counts for each position for this nt
  for (my $l=0;$l<=$#countsForThisNt; $l++) {#
    ## If the value for this position isn't defined,
    ## set it to 0.
    $countsForThisNt[$l]||=0;
    ## Print all the things
    printf "%.2f\t", $countsForThisNt[$l]/$., $l;
  }
  print "\n";
}

Save the script somewhere in your PATH, make it executable and run:
$ foo.pl file
A   0.13    0.40    0.13    0.00    0.00    0.40    0.47    0.07    0.40
C   0.40    0.27    0.00    0.00    0.00    0.07    0.07    0.13    0.27
G   0.20    0.20    0.67    1.00    0.00    0.47    0.20    0.73    0.20
T   0.27    0.13    0.20    0.00    1.00    0.07    0.27    0.07    0.13

Alternatively, if you're into the whole brevity thing, and enjoy some golfing, here's the same thing as a one-liner:
 perl -ne 'chomp;@F=split(//);for$i(0..$#F){$k{$F[$i]}[$i]++}}{for $nt(sort keys(%k)){print"$nt\t";for$i(0..$#{$k{$nt}}){$g=$k{$nt}[$i]||0; printf"%.2f\t",$g/$.}print"\n";}' file

Here an example in R using Biostrings and letterFrequency (as suggested by Devon Ryan).
library(Biostrings)

data <- read.table("DNA.txt", stringsAsFactors = F)

new <- matrix(nrow = 9, ncol = 15)

for(i in 1:9){
  for(j in 1:15){
    new[i,j] <- substring(data[j,], i, i)
  }
}

countTable <- matrix(nrow = 9, ncol = 4)
for(i in 1:9){
  columnSeq <- DNAStringSet(paste0(new[i,], collapse = ""))
  columnCounts <- letterFrequency(columnSeq, letters = "ACGT", OR = 0)
  countTable[i,] <- columnCounts
}
colnames(countTable) <- c("A", "C", "G", "T")

freqTable <- countTable/15
round(t(freqTable), digit = 2)
  [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
A 0.13 0.40 0.13    0    0 0.40 0.47 0.07 0.40
C 0.40 0.27 0.00    0    0 0.07 0.07 0.13 0.27
G 0.20 0.20 0.67    1    0 0.47 0.20 0.73 0.20
T 0.27 0.13 0.20    0    1 0.07 0.27 0.07 0.13

The whole new matrix generation can be written in a single line as do.call(rbind, strsplit(data, '')). The subsequent for loop can likewise be replaced.

Yes, there are many roads that lead to Rome.

benn’s answer works but is very un-R-like: you wouldn’t normally iteratively assign to items in a matrix: R has powerful vectorised operations which make this manual work unnecessary. Here’s a more idiomatic solution:
sequences = readLines('sequences.txt')
bases_matrix = do.call(rbind, strsplit(sequences, ''))

apply(bases_matrix, 2L, function (col) {
    str = DNAString(paste(col, collapse = ''))
    letterFrequency(str, letters = 'ACGT', OR = 0L, as.prob = TRUE)
})

This uses the same Bioconductor packages. Since this is such a simple problem, it can also be written without Bioconductor:
bases = strsplit(sequences, '')
# Use a data.frame here so we can use factors in the next step:
# R does not support matrices of factors. Ugh.
bases_by_column = setNames(do.call(rbind.data.frame, bases), seq_along(bases[[1L]]))
# Ensure that every column will be a complete set of ACGT frequencies
bases_by_column = lapply(bases_by_column, factor, c('A', 'C', 'G', 'T'))
sapply(lapply(bases_by_columns, table), prop.table)

Using modern R idioms from the ‘magrittr’ package, I’d write this as a pipeline; this very directly shows the sequence of transformations.
do.call(rbind.data.frame, bases) %>%
    setNames(seq_along(bases[[1L]])) %>%
    lapply(factor, c('A', 'C', 'G', 'T')) %>%
    lapply(table) %>%
    sapply(prop.table)

Double for loop is what I learned from java. In my opinion bioinformatics is about that it works, if you are only interested good code practices you should restrict to SO I guess.

@b.nota Nothing wrong with “it works”. But this being a Q&A site, I think we should definitely also teach best practices. in fact, there’s a problem in bioinformatics with this “it works” attitude, which leads to a proliferation of a lot of truly atrocious code. “But it works”. Yeah, just about … as long as you don’t nudge it.

I don't completely agree, no offense. With many bioinformaticians the BIO part is much more important than the best script practices. At the company where I work, they don't care about these best script practises, but do want someone that understands e.g., what a tissue resident memory cell is. That's just my point of view.

@b.nota It’s akin to lab (or generally science) best practices, to be honest: Many people use the same argument to justify sloppy lab work or dodgy use of statistics (p-hacking etc). At the moment bad software engineering is still somewhat accepted but I’m confident that this acceptance will lessen as the field matures. Just to clarify: this doesn’t so much concern your code, which works fine. But there’s a lot of genuinely bad, unmaintainable code in circulation in bioinfomatics tools that costs money and time to fix because of bad practices.

Yeah I agree with that part. I am PhD in biology, and learned to code with a java course here and a course in R there. I work with R almost daily, but unfortunately I still try to avoid apply functions (I do try them more often, but still the old school java double for loop is more intuitive for me).

Here is how you get the distributions by column using the TraMineR R package.
library(TraMineR)

sts.data <- c(
  "CAGGTAGCC",
  "CCGGTCAGA",
  "AGGGTTTGA",
  "TTGGTGAGG",
  "CAAGTATGA",
  "ACTGTATGC",
  "CTGGTAACC",
  "TATGTACTG",
  "GCTGTGAGA",
  "CAGGTGGGC",
  "TCAGTGAGA",
  "GGGGTGAGT",
  "TGGGTATGT",
  "GAGGTGAGA",
  "CAGGTGGAG"
)
seq <- seqdef(seqdecomp(sts.data,sep=''))

seqstatd(seq)

and here is the outcome of the seqstatd function
  [State frequencies]
   [1]  [2]  [3] [4] [5]   [6]   [7]   [8]  [9]
A 0.13 0.40 0.13   0   0 0.400 0.467 0.067 0.40
C 0.40 0.27 0.00   0   0 0.067 0.067 0.133 0.27
G 0.20 0.20 0.67   1   0 0.467 0.200 0.733 0.20
T 0.27 0.13 0.20   0   1 0.067 0.267 0.067 0.13

   [Valid states]
  [1] [2] [3] [4] [5] [6] [7] [8] [9]
N  15  15  15  15  15  15  15  15  15

   [Entropy index]
   [1]  [2]  [3] [4] [5]  [6]  [7]  [8]  [9]
H 0.94 0.94 0.62   0   0 0.78 0.87 0.62 0.94


text:
Pipeline (Pypiper) returns error, but individual command works fine
I am running a script file on a cluster with slurm through sbatch. Locally it runs perfectly well on 2 different linux versions (linux mint, ubuntu 12.04), however, it returns -9 error code while running on a cluster which is CentOS 7.3. When we are using Pypiper and run the whole .py script file, it fails when it starts running trimmomatic program
http://www.usadellab.org/cms/?page=trimmomatic
However, when running just trimmomatic program step by itself it does not fail and finishes successfully. We have a suspicion that something may be wrong with how Pypiper runs on the cluster, however, feeling rather stuck.
Any suggestions would be greatly appreciated.
Here is the error that is printed out:
ILLUMINACLIP: Using 2 prefix pairs, 22 forward/reverse sequences, 0
forward only sequences, 0 reverse only sequences
</pre>
Process 109391 returned: (-9). Elapsed: 0:00:51. Peak memory: (Process:
5.211GB; Pipeline: 0.0GB)

### Pipeline failed at:  (12-14 15:50:47) elapsed:0:23:09 _TIME_

Total time: 0:23:09

Change status from running to failed
Traceback (most recent call last):
  File "...../salmon_rna_seq.py", line 344, in <module>
sys.exit(main())
File "...../salmon_rna_seq.py", line 304, in main
pipeline.process()
File "...../salmon_rna_seq.py", line 277, in process
self.trim_reads()
File "...../salmon_rna_seq.py", line 195, in trim_reads
self.pipe_manager.run(cmd, self.sample_files['r1']['fastq']['trimmed_paired'])
File "...../pypiper.py", line 508, in run
process_return_code, local_maxmem = self.callprint(cmd, shell, nofail, container)  # Run command
File "...../pypiper.py", line 691, in callprint
self.fail_pipeline(e)
File "...../pypiper.py", line 1111, in fail_pipeline
raise e
OSError: Subprocess returned nonzero result.

Pypiper terminating spawned child process 108475...
child process terminated
slurmstepd: error: Exceeded step memory limit at some point.

The last line may be informative. Although, we are still not sure what we could possibly do.

Slurm is killing the process because it's using more memory per core than the default allocation (that's what slurmstepd: error: Exceeded step memory limit at some point means). I've never used pyPiper (like seemingly most people, I use snakemake), but presuming that you have it running either srun or sbatch at some point then specify something like --mem-per-cpu 10G.

Actually I tried that and it still does not work. Will try to write to the cluster support

You might need a larger value than 10 gigs. You can use scontrol to see the exact threshold getting crossed that's causing slurm to kill the job.

Awesome. We increased it to 20Gb and it run successfully.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/

sort: votes


text:
changing color key range to specified range in heatmap.2 function
I have a tab separated text file as shown below
a   0.311   0.510   0.123   0.002   0.001   0.417   0.572   0.074   0.169
b   0.324   0.592   0.070   0.028   0.028   0.535   0.535   0.127   0.113

I am trying to use heatmap.2 function from gplots package and want to change the color key so that it ranges from 0 to 1
Following is my code:
d <- read.table("file.txt",sep="\t",header=F,row.names = 1)
colnames(d) <- c("-3","-2","-1","+1","+2","+3","+4","+5","+6")
d_matrix <- data.matrix(d)
colors<-seq(0,1,length.out=100)
heatmap.2(d_matrix,scale="none",col = redblue(75),trace="none",Colv = FALSE,Rowv = FALSE,breaks=colors)

But this gives me error:
Error in image.default(1:nc, 1:nr, x, xlim = 0.5 + c(0, nc), ylim = 0.5 +  :
  must have one more break than colour

I want to change to color key such that it ranges from 0 to 1. How can I do it in the function heatmap.2
I added breaks argument to the heatmap.2 function as suggested from this post here but no luck..

You just need to add more colors. The following, for example, will work:
heatmap.2(d_matrix, scale="none", trace="none", Rowv=F, Colv=F, breaks=seq(0,1,0.01), dendrogram="none", col=redblue(100))

Note that length(seq(0, 1, 0.01)) is 101, therefore redblue() must return one fewer (100) colors (or colours if you prefer extra vowels).

This worked. I have 2 more issues. I am having some sky blue lines in the color key which I don't want and secondly my rownames a and b are very big. Is there a way to reduce the size?? Thanks

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
how do I predict bacterial small non-coding RNA for a specific mRNA?
I am working on Vibrio parahaemolyticus. I have a gene that it is possible might be regulated by small non-coding RNAs. How can I predict possible sRNAs that target the transcript of this gene. Most software offers the opposite of what I want, i.e., they predict targets for small RNAs!
I'm trying to predict bacterial small non-coding RNA in general (Support Vector Machine) but I got more than 8 thousand sequences within intergenic regions alone. Is there a code for aligning these sequences to my transcript?

Answer from @llrs, converted from comment:
You could try to download RNAseq data from this bacteria and see if some of the RNA align to your sequence. You could try to blast your gene and look for the sequences that have low identity similarity (or classify them by type and then look what is each hit).

Good idea but from my experience, normal alignment is not optimized for sRNA. Tools like targetRNA2 are optimized for sRNA alignment such as imperfect pairing and energy ..etc.

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Using Python, how to convert a pandas DataFrame into a VCF?
Let's say I have a pandas dataframe with fields CHROM, POS, ALT, REF. In this special case, I also wouldn't care about IDor FILTER, INFO could be blank (or meaningless) and we'll write QUAL as each 40.
import pandas as pd
example_dict = {'CHROM':[20, 20, 20, 20, 20], 'POS':[14370, 17330, 1110696, 1230237, 1234567], 'REF':['G', 'T', 'A', 'T', 'GTCT'], 'ALT':['A', 'A', 'G', 'C', 'G']}

df = pd.DataFrame(example_dict)
print(df)

   CHROM      POS   REF ALT
0     20    14370     G   A
1     20    17330     T   A
2     20  1110696     A   G
3     20  1230237     T   C
4     20  1234567  GTCT   G

In this case, saving the output with the extension .vcf via df.to_csv("myfile.vcf") with the other required fields would created the tab-delimited format of the VCF, but it wouldn't include the mandatory header. https://samtools.github.io/hts-specs/VCFv4.1.pdf
##fileformat=VCFv4.1
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=file:///seq/references/
...
#CHROM POS ID REF ALT QUAL FILTER INFO

The most popular package to work with VCF's appears to be PyVCF: https://pypi.python.org/pypi/PyVCF
However, one needs to begin with a VCF as a template, and add to this (from the documentation):
vcf_reader = vcf.Reader(filename='vcf/test/tb.vcf.gz')
vcf_writer = vcf.Writer(open('/dev/null', 'w'), vcf_reader)
for record in vcf_reader:
     vcf_writer.write_record(record)

It appears the only information necessary from vcf.Reader(filename='vcf/test/tb.vcf.gz') is the metadata?
print(vcf_reader.metadata)
OrderedDict([('fileformat', 'VCFv4.0'), ('fileDate', '20090805'), ('source', ['myImputationProgramV3.1']), ('reference', '1000GenomesPilot-NCBI36'), ('phasing', ['partial'])])

Otherwise, tb.vcf.gz already has variants, which I think means I would need to feed in a VCF with no variants?:
for record in vcf_reader:
     print(record)

Record(CHROM=20, POS=14370, REF=G, ALT=[A])
Record(CHROM=20, POS=17330, REF=T, ALT=[A])
Record(CHROM=20, POS=1110696, REF=A, ALT=[G, T])
Record(CHROM=20, POS=1230237, REF=T, ALT=[None])
Record(CHROM=20, POS=1234567, REF=GTCT, ALT=[G, GTACT])

You will have to define the header either from an existing VCF file or hardcoded into you python script.  Then write the header the the output VCF file then write the dataframe to the same file with the mode options set to 'a' to append to the end of the file.
You will need to order the dataframe columns to match the VCF spec and insert the empty/ default ID, FILTER, INFO and QUAL columns
import pandas as pd

example_dict = {'CHROM':[20, 20, 20, 20, 20], 'POS':[14370, 17330, 1110696, 1230237, 1234567], 'REF':['G', 'T', 'A', 'T', 'GTCT'], 'ALT':['A', 'A', 'G', 'C', 'G']}

df = pd.DataFrame(example_dict)

df['ID'] = ""
df['QUAL'] = 40
df['FILTER'] = ""
df['INFO'] = ""

df = df[['CHROM', 'POS', 'ID', 'REF', 'ALT', 'QUAL', 'FILTER', 'INFO']]

header = """##fileformat=VCFv4.1
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=file:///seq/references/
#CHROM POS ID REF ALT QUAL FILTER INFO
"""

output_VCF = "myfile.vcf"
with open(output_VCF, 'w') as vcf:
    vcf.write(header)

df.to_csv(output_VCF, sep="\t", mode='a', index=False)

Then myfile.vcf looks like
##fileformat=VCFv4.1
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=file:///seq/references/
#CHROM POS ID REF ALT QUAL FILTER INFO
CHROM   POS ID  REF ALT QUAL    FILTER  INFO
20  14370       G   A   40
20  17330       T   A   40
20  1110696     A   G   40
20  1230237     T   C   40
20  1234567     GTCT    G   40

Thanks for the help. BTW, is fileDate year-month-day or year-day-month? I'm guessing there's no standard.

fileDate is not a required header field you can remove it it you want


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

Locus1:
>Indiv1
ACGT
>Indiv2
ACGC
>Indiv3
ACGA

Locus2:
>Indiv1
GACT
>Indiv2
GACA
>Indiv3
GAGT

Locus3:
>Indiv1
GGTC
>Indiv2
GGCT
>Indiv3
GCTA

(These example alignments leave out quality scores for simplicity, assume that each has high quality across the board.)
What algorithm or system would be used to graph them on a PCA? How many dimensions would that PCA have? Is it possible to model this by hand so that I can understand what's going on/how it is built?

Can you update your post to somewhere mention that this is related to population structure? There are other uses of PCA using BAM files that use completely unrelated methods (e.g., deepTools).

In general, the process works as such:

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

@terdon "...according to the number of copies of the reference allele that it has (so, 0-2)."    So the PCA is related to heterozygosity?

You said each locus was scored often based on number of copies- what else could it be scored on?

So you build a matrix (often) based on allele copy (0, 1, or 2) then make a covariance matrix and do associated eigendecomposition stuff to generate a PCA?

One could consider making a numeric score for the genotype, such that each non-reference genotype got a different value (instead of them being all lumped into 0/1). I've not seen someone try that, so I don't know whether it's useful. Usually PCA is done via SVD these days, but really however you want to do it should lead to similarish results.


text:
Convert rs ID of one hg build to rs IDs of another build
I have a list of dbSNP rsIDs for GRCh37 and I want to convert them to the equivalent IDs in GRCh38. This is using the most recent dbSNP build (150 as of the time of this post). Is there any ID mapping available? If not, are there any tools I can use?

Did you try the liftOver tool of UCSC? You'll need a BED file with your SNPs coordinates for this.
You can also do this in R with rtracklayer.
library(rtracklayer)
?liftOver

There's no need to change the coordinates, the question is regarding the IDs themselves.

What I have is a  list of dbSNP IDs. I have downloaded source("http://bioconductor.org/workflows.R") workflowInstall("liftOver"), but not sure if ths will accept just a list of IDs.

No it won't accept only rsIDs, so you'll need to get the coordinates first with e.g., biomaRt. Sorry I understood that you had coordinates of the SNPs.

Thanks @b.nota I was able to get coordinates from biomaRT.

You can assume that the overwhelming majority of rsIDs are the same between GRCh37 and GRCh38 (they're semi-stable IDs). There are, however, a number of rsIDs that are present only in GRCh37, which you can find here. Note that the format of this file is a bit strange, it's chromosome|position|ID|weight, where position is sometimes empty, weight is typically 1, and ID is the rsID without the rs prefix.
Also, some rsIDs have a different strand. This is less likely to be an issue for you, but it's there in case you need it.

Thanks so much @Devon Ryan for all the useful information.

See also CrossMap:
http://crossmap.sourceforge.net/
20, 21, 22, 23, 24, 30

These are only the base positions that have changed between GRCh37 and GRCh38, not the IDs. For example if a study (e.g. found via pubmed) refers to a particular rsID, the genome build does not matter.
rsIDs do not depend on a reference genome. They point to a specified locus regardless of the differences in genomic assemblies. As a result, you cannot convert rsIDs of one genome build to rsIDs of another build. And there is no such mapping tool available, even if they pretend to be. rsIDs are assigned by dbSNP maintainers, and are stable. The numbers may only change on merge, where multiple numbers were assigned at the same genomic location, and the higher number gets merged into a lower number by the dbSNP maintainers, unless there is a very well-published, notable rsID mentioned in the press, in this case that merged number is considered the mentioned one, not the smaller one.

RefSNP number is the stable accession regardless of the differences in genomic assemblies.

For more information go to the “RefSNP Number Stability” section at https://www.ncbi.nlm.nih.gov/books/NBK21088/
The rsIDs may also be deleted in some cases see the following links on the explanation:

https://www.ncbi.nlm.nih.gov/books/NBK44467/#Build.i_found_rs123456_in_the_current_bu
https://www.ncbi.nlm.nih.gov/books/NBK44468/#Build.the__rscurrent_field_in_rsmergearc

Thus, specifying the location by rsID is stable, and can be only affected by an eventual merge or delete from the dbSNP database, if any, not by a new genome build. You can only “update” the rsIDs from time to time to make sure they all are the main rsID, not those larger numbers already merged to a smaller one. It should be noted that the base positions (offsets) of the same genomic location may change between the builds of the reference genome. And they did surely change between GRCh37 and GRCh38.
If your list is relatively small, you can get most recent rsIDs for your chromosome and offset via dbSNP API. It is available at https://www.ncbi.nlm.nih.gov/snp/docs/eutils_help/
The query for the API is the same as for the “advanced search” of dbSNP at https://www.ncbi.nlm.nih.gov/snp/advanced
If you specify data in a “Base Position” line, it will look for GRCh38 base position. If you need to specify the base position for a GRCh37, use “Base Position Previous” line. You may also need to add NOT MergedRS[All Fields] to exclude rsIDs that have already been merged.
For example, to get the rsID for a chromosome 1 GRCh37 base position 249222527, the API query URL is the following https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=snp&term=1[Chromosome]+AND+249222527[Base+Position+Previous]+NOT+MergedRS[All+Fields]
You can do queries from a script if the number of requests is relatively small. You may find a perl script, for example, that does the job.
Even though the rsID numbers can be merged or deleted, and the changes are reflected in a new dbSNP database build, given the fact that the changes in numbers do not depend on the genome build, is better to always rely to the rsID number from the latest dbSNP database build.
The changes in rsID numbers are changed due to merges or deletions between dbSNP builds relatively infrequently, given the total number of the IDs in the database.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Filter out PCA outliers automatically
I am new to bioinformatics and PCA. What I am trying to do is to remove bad cells from a dataset that was obtained with scRNA-seq for 46 000 of genes (rows) vs 68 000 of cells that are of different types of C. Elegans that was generated in the following paper.
After running PCA:
fontsize <- theme(axis.text=element_text(size=12),
axis.title=element_text(size=16))
plotPCA(sce, pca_data_input="pdata") + fontsize

I see the following chart:

I suppose that the cells that are going to the left bottom corner are the bad ones since they do not fall into some restricted region. Are there any R program that can automatically remove the bad cells leaving me with just the majority of cells clustered within an ellipse? I know there is identify method that allows for clicking and labelling manually cells that I think may be outliers (see this) but I do not like it. I would prefer some program that would do that for me, giving some statistical output of the reliability of the removal of the cells.

I suspect something has gone wrong if 96% of your gene expression variance is in the first two components (and 87% in the first component). Have you normalised counts for each cell? Do you have the gene expression matrix oriented correctly? Have you transformed the expression to be in log space? Are you scaling the PCA such that each gene contributes the same amount to the variance?

Thank you very much for the advice! I eliminated low-count genes, normalized the dataset, eliminated the low-count cells and PCA started to look much better in terms of PC1 and PC2 percentages, although still far from perfect and there is no separation into unique groups: https://bioinformatics.stackexchange.com/questions/3232/improve-scrna-seq-dataset-for-further-analysis

I still do not know the answer to this question

It is not an easy task to automatically identify outliers in a multivariate space. Have you tried to look at this or this ?


text:
Filter out outliers of the scRNA-seq (heterogenous cells)
I am new to data science. I have a dataset of single-cell gene expression from multiple cell types in C. Elegans. The dataset is from the paper Comprehensive single-cell transcriptional profiling of a multicellular organism
My main question is, Which approaches should I use for filtering out bad cells in this case when we have multiple cell types in the dataset?
So far I tried to filter out genes that have too high mitochondrial genes content following the Bioconductor “simpleSingleCell” workflow.
However, the tutorial specifically says that the method of filtering out based on mitochondrial genes most probably will not work when the dataset has multiple cell types:

Analyzing all cell types together would unnecessarily inflate the MAD and
      compromise the removal of low-quality cells, at best; or lead to the
      entire loss of one cell type, at worst.

Any suggestions would be greatly appreciated.

By filtering out what do you mean: removing the sample from the analysis when there are too many mithocondrial genes or removing the mitochondrial genes ? (This seems important in scRNAseq and I haven't worked with it, so excuse me if it is a naïve question)

removing the bad quality samples

Not a standard, but I found this course helpful http://hemberg-lab.github.io/scRNA.seq.course/ . I personally filter cells with too low or too high transcript counts, cells with low count of detected genes, and cells with high spike-in/endogenous RNA ratio. Prior to this pipeline I also filter out ribosomal RNA. Take it with a grain of salt because I'm still experimenting on it. Also, this answer may help https://bioinformatics.stackexchange.com/a/3171/1771

From what I known, there is no clear consensus in the field and it depends on the type of cells you are interrogating.
However, if mitoRNA/endogenousRNA ratio does not fit for your purposes, other option is to check at the total number of genes/transcripts detected in each cell. In this way, you can filter out cells that have considerable less genes/transcripts detected than the rest since this may be indicative of bad quality cells for any reason (apoptosis, RNA degradation, sequencing itself, etc). For instance, to put a threshold of >1000 genes/transcripts in a cell to be considered for further analysis.
Consider also to take a look at the number of mapped reads for each cell, as cells with low mapped reads respect to the rest could be potentially problematic.

What does it mean ‘mapped reads’? How is it different from the expression matrix itself?

If you are using UMIs, usually I filter out those that are supported by less than 3 reads and some that do not fulfill QC requirements. In this case UMIs and mapped reads does not correspond (although should be very similar). At the end is another way of perform QC of the data but at an earlier stage of the analysis (before counting molecules and assign them to genes).


text:
Normalizing RNAseq for PCA and CCA
Usually the expression data is transformed to log space using either RPKM, FPKM or CPM, this is required when looking for differential expression because the data is tested against the normal distribution(limma) or the negative bionimal distribution (DESeq2).
In weighted gene co-expression analysis (WGCNA) there is also the recommendation to normalize the data (see FAQ 4) but says "RPKM, FPKM, or simply normalized counts doesn't make a whole lot of difference for WGCNA" because it is based in correlations.
For PCAs the data should be on the same scale to avoid missing the scale as a factor. See this for more references. But should it be normalized too to look like a normal distribution using logarithms?
How does this extend to other methods like CCA? Should the data look like a normal distribution for CCA? As it is based on correlation I don't expect much changes, so I am using raw counts (RSEM output), but the examples I've seen use normalized data.
This question comes after reading a comment

I've given an answer about using how I normalise for PCA here. PCA works best if the data are normal-ish, but as models go it is fairly robust to fairly large deviations from normal. Transforming is not essential, but if data can be easily transformed into something that looks a bit more normal, then conclusions derived from looking at the PCA should be a bit more robust.
I'm not familiar with canonical-correlation analysis (CCA), but because I see sums and square roots in the formulas, I would assume it has similar data cleaning requirements to PCA. Unfortunately the R documentation for cancor doesn't have any information about any assumptions made on the data.


text:
Finding a single open reading frame with ribosomal binding site, using Biopython
I'm given a Fasta file, containing a large DNA(over 115,000 long) sequence, and I am tasked with finding a single large open reading frame contained within the DNA sequence using Biopython.
I'm aware this has been asked before consistently however i can't find anything for actually finding the source of my ORF, but apologies for any overlap.
Now from other sources and the Biopython cookbook I've translated my sequence and found six open reading frames (three for each strand) and their positions within the sequence;
def find_orfs_with_trans(seq, trans_table, min_protein_length):
    answer = []
    seq_len = len(seq)
    for strand, nuc in [(+1, seq), (-1, seq.reverse_complement())]:
        for frame in range(3):
            trans = str(nuc[frame:].translate(trans_table))
            trans_len = len(trans)
            aa_start = 0
            aa_end = 0
            while aa_start < trans_len:
                aa_end = trans.find("*", aa_start)
                if aa_end == -1:
                    aa_end = trans_len
                if aa_end-aa_start >= min_protein_length:
                    if strand == 1:
                        start = frame+aa_start*3
                        end = min(seq_len,frame+aa_end*3+3)
                    else:
                        start = seq_len-frame-aa_end*3-3
                        end = seq_len-frame-aa_start*3
                    answer.append((start, end, strand,
                                   trans[aa_start:aa_end]))
                aa_start = aa_end+1
    answer.sort()
    return answer

orf_list = find_orfs_with_trans(record.seq, table, min_pro_len)

for start, end, strand, pro in orf_list:
print("%s...%s - length %i, strand %i, %i:%i" \
      % (pro[:30], pro[-3:], len(pro), strand, start, end))

OUTPUT
TNRQVYGGTLQSLRTGTGIYSRLASSPTNR...CEA - length 160, strand 1, 7950:8433
 GIPPGRTEGLGRYVHGESIFLEATLVPEPQ...VQS - length 171, strand -1, 19275:19791
 ISVHLRRYFSVLLRAPVALNADSRVAGPLD...ITP - length 190, strand 1, 34079:34652
 GCVNGNFPDYRVAVDDPGALVVGGELGQAL...VNT - length 771, strand -1, 39335:41651
 LLKLLASRLTWLLVLSAALGFLTSVSYRLG...SPG - length 235, strand -1, 39358:40066
 PGDVIVSKVPADNLRPRMSEINDFTNSIII...VLV - length 826, strand 1, 39362:41843
 GPGRQSSSAYERDQRLHQQHHHQRERHQVR...QRS - length 764, strand 1, 39385:41680
 LLLSSTSSASCWLKKFFFRLLLLSRSVALL...LAL - length 158, strand -1, 40459:40936
 TRLIIRPPSHTTLYYDPYPQPLKVSSSVPC...RIT - length 159, strand -1, 54663:55143
 IGAFSRGYNQGWVLIEAQEVRTTSPLFPRS...RDQ - length 201, strand 1, 56104:56710
 KARGSSLGAIYCRRLKRRSLPSWQIESFNC...TNS - length 150, strand -1, 62497:62950
 RRPARAGAELNKDDIRDTCLLSCSEVDNKV...VTL - length 168, strand 1, 72339:72846
 SGPTTVRTSAPQRCWRTSLKPKLSLSDSTE...PVA - length 177, strand -1, 114494:115028
 ERELWLETGPPTPLWGTGPDCSRAALNRQR...PRR - length 157, strand 1, 114953:115427
However I'm unsure if these open reading frames contain the ribosome binding site(consensus sequence of AAGGAGGTG and occurring between 6 to 9 nucleotides upstream from the positive 5' strand). Is there anyway to add to my code to include this? or would it already be included with Biopython?
I'm left with the output however I've no idea how to link this with possibly an online source like NCBI to actually find the source of my ORF (i.e. which protein it comes from). Is there anyway to use Biopython for this as well?


text:
Show presence of known mutation in RNA-seq data
We have RNA-seq fastq data from control (WT) patients and a patient with a point mutation at a known location in one gene.
I'd like to retrieve the reads aligning to that gene and show the presence of the mutation.
I can think of 2 strategies:
A)
Alignment and visualization, based on GATK:

Run STAR 2-pass
Extract alignments of interest, e.g.:
samtools view input.bam "Chr10:18000-45500" > output.sam
Visualize:

IGV
IGB
ASCIIGenome (see also this SE answer)

How can I count the reads showing mutation vs wild type?

B)
The second approach would simply look for the count numbers mapping (e.g. salmon) to the transcript ID(s) corresponding to mutation and WT, based on this SE answer.
Is there a better way to do this?

For counting reads I use mpileup, e.g. samtools mpileup --reference hg38.fa -r Chr10:18000-45500 input.bam, which will give base-resolution coverage for a BAM file.
I've written my own script to process mpileup output and make it easier to understand. By default it reports read coverage as a proportion of total coverage, but this can be modified by using the -counts command line argument:
$ samtools mpileup -r tig00018708_tig00000379:210665-240664:2107-2130 --reference trimMmerged_tig00018708_other.fasta local_Sampled_50M_vs_trimMmerged_tig00018708_other.bam  | /bioinf/scripts/readstomper.pl -counts
[mpileup] 1 samples in 1 input files
<mpileup> Set max per-file depth to 8000
Assembly,Position,Coverage,ref,cR,pR,A,C,G,T,d,i,InsMode
tig00018708_tig00000379:210665-240664,2107,10,C,10,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2108,11,G,11,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2109,11,T,11,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2110,11,T,11,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2111,12,G,12,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2112,14,G,14,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2113,13,C,13,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2114,13,G,13,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2115,13,C,2,0.15,0,0,0,11,0,0
tig00018708_tig00000379:210665-240664,2116,15,T,14,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2117,15,G,14,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2118,18,A,17,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2119,19,T,18,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2120,19,C,18,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2121,18,A,7,0.41,0,0,10,0,0,0
tig00018708_tig00000379:210665-240664,2122,18,G,18,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2123,18,T,18,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2124,18,C,18,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2125,18,C,9,0.50,0,0,0,9,0,0
tig00018708_tig00000379:210665-240664,2126,18,C,14,0.78,2,0,1,1,0,0
tig00018708_tig00000379:210665-240664,2127,18,C,18,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2128,18,T,15,0.83,0,0,3,0,0,0
tig00018708_tig00000379:210665-240664,2129,17,A,13,0.76,0,0,4,0,0,0
tig00018708_tig00000379:210665-240664,2130,16,A,9,0.56,0,6,1,0,0,0

Note that only non-reference coverage is reported in the A/C/G/T/d/i columns. I have found that this makes it easier for subsequent data analysis, but your mileage may vary.

If you only have one gene and you only need to do this once then the simplest possible workflow is to generate the aliment using STAR (optimally with the two pass method) and open the two resultant .bam files in IGV, the coverage column at the stop above your alignment track should clearly show the counts of reference and non-reference supporting reads.  The IGV documentation shows a screen shot the pop-up coverage indicator and allele specific read count.

SCmut is a method to detect cell-level mutation from single-cell RNA-sequencing, published by Vu et al. (2019)
https://github.com/nghiavtr/SCmut


text:
Large dataset normalization for PCA
I need to normalize a large (400Mb) dataset for doing PCA analysis. I want to use scran for doing that:
biocLite("SingleCellExperiment")
biocLite("scran")
library(SingleCellExperiment)
library(scran)
list_of_sce <- list()
# Looping though the UMI_count 'split_factor' columns at a time
split_factor = 500
for(i in seq(1,ncols, split_factor)) {
  num_loop = floor(i / split_factor) + 1
  idx = ncols
  if (i + split_factor < ncols) {
  idx = i + split_factor
  }
  sce <- SingleCellExperiment(list(counts=UMI_count[, i : idx]))
  # Normalization of the dataset containing
  # heterogenous cell data (different cell types)
  clusters <- quickCluster(sce) ## <- ** error appears here **
  sce <- computeSumFactors(sce, cluster=clusters)
  list_of_sce[[num_loop]] <- sce
}
cbind(list_of_sce)

However, when I try to run this it produces an error at the quickCluster step:
Error: cannot allocate vector of size 23.7 Gb

It feels like cbind at the end is not the right way to go, and I need to somehow combine sce and re-normalize them again somehow. The question is how to do that in a correct way: normalizing large dataset for PCA analysis when it can not be fully loaded into memory?

BTW, there are many datasets on kaggle.com

Which error do you get when you normalize the whole dataset? How big is your dataset?

I normalized a high throughput dataset for a school project using DESeq library using the script bellow. The code is based on a lesson I had. My goal was determine the over expressed genes, but the normalization step should be the same.
#####################################################################
### DATA PRE-PROCESSING: Normalising gene expression distributions
#####################################################################
# Genes that are not expressed at a biologically meaningful level
# in any condition should be discarded to reduce the subset of genes
# to those that are of interest, and to reduce the number of tests carried out
# downstream when looking at differential expression.
# Using a nominal CPM value of 1 (which is equivalent to a log-CPM value of 0
# a gene is deemed to be expressed in a given sample if its transformed count
# is above this threshold, and unexpressed otherwise.
# Genes must be expressed in at least one group (or roughly any three samples
# across the entire experiment) to be kept for downstream analysis.

library("DESeq")

keep.exprs <- rowSums(cpm>1)>=3 # is filtering is ok
DGEfiltered <- DGEobject[keep.exprs,, keep.lib.sizes = FALSE]

## Normalisation generally required to ensure that expression distributions
## of each sample are similar across the entire experiment.
## Any plot showing the per sample expression distributions, such as a density or boxplot,
## is useful in determining whether any samples are dissimilar to others.
## Normalisation by the method of trimmed mean of M-values (TMM)
## is performed using the calcNormFactors function in edgeR.
# cpmF = filtered subset # log.cpmF = log cpm filtered subset
# Copy unnormalized, filtered data to show in boxplot
DGEf.notnorm <- DGEfiltered # normaliseerde
DGEf.norm <- calcNormFactors(DGEfiltered, method = "TMM")
DGEf.norm$samples$norm.factors

# --> scaling factors relatively close to 1 --> mild TMM normalisation

## Boxplots of log-CPM values showing expression distributions
## for unnormalised data (A) and normalised data (B) for each sample
par(mfrow = c(1,1))
par(mar = c(7,4,2,1))
# Boxplot unnormalised data
DGEf.notnorm$samples$norm.factors
# [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1
log.cpmF.notnorm <- cpm(DGEf.notnorm, log = TRUE)
boxplot(log.cpmF.notnorm, las = 2, col = sample.col, cex = 0.9,
        main = "A. Unnormalised data", ylab = "Log-cpm")
# Boxplot normalised data
DGEf.norm$samples$norm.factors

log.cpmF.norm <- cpm(DGEf.norm, log = TRUE)
boxplot(log.cpmF.norm, las = 2, col = sample.col, cex = 0.9,
        main = "B. Normalised data", ylab = "Log-cpm")
# --> only when scaling factors are not close to 1 you will see a big difference

Hello nemo, thanks for your answer, looks good. Could you add also some info where the script comes from? Maybe a reference or at least a bit of reasoning? I also believe that you should include some libraries.

I wrote it for the most part my self for a school project. The rest of the code is from the lesson I had. It was for an high throughput analysis on dataset for a disease. To determine the over expressed genes.

Cool, if the materials of the class are available it would be good to link it, but that sounds already quite good to me. I tried to [edit] your post (& add there the info), check if it still makes sense and do not hesitate to make a further edit to fix it or to add more details/links.


text:
Error in .computeSumFactors: cells should have non-zero library sizes
While trying to run the following code in R:
list_of_sce <- list()
# Looping though the UMI_count 'split_factor' columns at a time
split_factor = 500
for(i in seq(1,ncols, split_factor)) {
  num_loop = floor(i / split_factor) + 1
  idx = ncols
  if (i + split_factor < ncols) {
    idx = i + split_factor
  }
sce <- SingleCellExperiment(list(counts=UMI_count[, i : idx]))
# Normalization of the dataset containing
# heterogenous cell data (different cell types)
clusters <- quickCluster(sce)
sce <- computeSumFactors(sce, cluster=clusters)
list_of_sce[[num_loop]] <- sce
}

I am getting the error:
Error in .computeSumFactors(assay(x, i = assay.type), subset.row =
subset.row,  :
cells should have non-zero library sizes

There are warnings, so when I run warnings() command I see the following:
> warnings()
Warning messages:
1: In .computeSumFactors(assay(x, i = assay.type), subset.row =
subset.row,  ... :
encountered negative size factor estimates
2: In .computeSumFactors(assay(x, i = assay.type), subset.row =
subset.row,  ... :
encountered negative size factor estimates
3:...
...
39:....

I do not understand why is this happening, and how to debug it. Any suggestions would be greatly appreciated.

What is the library size of the samples? computeSumFactors is failing because it can't compute the library size. I am not sure if your code is working as intended, shouldn't the sce be inside your loop?

Do you have a sample with all 0s?

Your code is different from another question you posted (here). Is this the same data/problem as in the linked question?

It is the same. Where is the difference? Can't see it. These are related but different questions, so I asked them separately.


text:
How to compute the chance of failing to detect a gene given the detection limit of a protocol
In Shapiro et al., when discussing about loss of molecules as source of error in single-cell sequencing, it is written that:

Another source of error is losses, which can be severe. The detection limit of published protocols is $5$–$10$ molecules of mRNA. If, as seems likely, the limit of detection is primarily determined by losses during sample preparation, this would indicate that $80$–$90\%$ of mRNA was lost. Or, to put it the other way around, a $90\%$ loss leads to an approximately $50\%$ chance of failing to detect a gene that is expressed at a level of seven mRNA molecules (from the binomial distribution).

How is this probability computed using the binomial distribution? I thought that $90\%$ loss corresponds to $5$ detected molecules, and I assume that $k=7$ for the binomial calculation, but I am unable to go further.

A 90% loss can be rephrased as a 10% chance of detecting anything. So what we want to find is the probability of detecting 0 molecules, when we start with 7 and have 10% probability of success. Once can do that in R as follows:
> pbinom(0, 7, 0.1)
0.4782969

So ~50%, as they stated. I suspect that part of the confusion arises from the fact that the detection limit is due to lowly-expressed genes/transcript being heavily affected by this loss. So the probability of detecting a single molecule out of 4 original molecules (assuming 90% mRNA loss) is 34%, for 3 molecules it's 27%, for 2 it's 19% and for 1 it's 10%. I think the threshold of 5 is mentioned more because it's a nice round number than there's anything particularly different in detecting a gene with 4 vs. 5 molecules (34 vs. 41% probability).

I was reading again your answer to understand the first sentence. If I understood correctly, the mapping from $90%$ loss to $10%$ chance happens computing the probability of detecting a single molecule out of $1$ original molecule. Right?

Correct, since there's a 90% loss of all signal, there's a 10% probability of detecting any given molecule.

Detecting zero molecules out of 7 can be rephrased as failing to detect a molecule, 7 times independently. The probability of this is 0.9 to the power of 7. Python tells me (0.9 ** 7) that this is 0.4782969. That's probably what the modelization in terms of binomial distribution means, but to me it is much easier to understand in terms of this very basic probability calculation.


text:
How to indicate the END of a haplotype block in VCF?
In VCF I know how to indicate that two genotypes are in the same phase by using consecutive "0|1" and "1|0" genotype fields, for example. However, how do I deal with the case that the first two genotypes are in phase and the second 2 are in phase, but there is no phase between the 2nd and 3rd?
According to the format I think I need to write:
chr1 100 . A T . . . GT 1|0
chr1 200 . A T . . . GT 0|1
chr1 300 . A T . . . GT 1|0
chr1 400 . A T . . . GT 0|1

but that implies that the SNP at 200 and the SNP at 300 are in phase, when in fact they are not.
How to indicate the END of a haplotype block in VCF?

According to the VCF specification, haplotype blocks are defined by the "Phase Set" [PS] tag:

PS : phase set. A phase set is defined as a set of phased genotypes to which this genotype belongs. Phased genotypes for an individual that are on the same chromosome and have the same PS value are in the same phased set. A phase set specifies multi-marker haplotypes for the phased genotypes in the set. All phased genotypes that do not contain a PS subfield are assumed to belong to the same phased set. If the genotype in the GT field is unphased, the corresponding PS field is ignored. The recommended convention is to use the position of the first variant in the set as the PS identifier (although this is not required). (Non-negative 32-bit Integer)

A good example of this can be found on a 10X Genomics support page about phasing:
CHROM   Pos REF ALT GT  PS
chr1    1000    A   C   0|1 1000
chr1    1010    T   G   1|0 1000
chr1    2000    C   T   0|1 2000
chr1    2005    T   G   0/1 2000
chr1    2008    G   C   0|1 2000

In this example we have two phase blocks, denoted by PS=1000 and PS=2000. PS=1000 spans position 1000-1010, and PS=2000 spans position 2000-2008. In PS=1000, haplotype 1 contains the REF A allele at position 1000, and the ALT G allele as position 1010, while haplotype 2 contains the ALT C allele at position 1000 and the REF T allele at position 1010.

Note that this phase set identifier needs to be a single numeric value, which means that sub-blocks and overlapping blocks cannot be defined by this field.

Two options that immediately come to mind, depending on what you hope to accomplish down the pipeline: you might simply mark the unphased sites with a "/" instead of "|" to indicate they are not phased; or, if you have multiple haplotype blocks that are all independent of each other, you could use the chromosome field to identify each unique haplotype (e.g. chr1h1, chr1h2, chr2h1, chr2h2, etc).

There is a new idea of reporting the "phased haplotype/genotype block" with unique block ID.
phaser uses "PI" (unique phase-index) and "PG" (phased-genotype).
GATK is also using this idea of "unique-block" and "phased-genotype"; See the "HP" tag.
Personally, I feel the representation of "phaser" is way better and comprehensive.
I would represent your data as:
chr1 100 . A T . . . GT:PI:PG    1/0:3:1|0
chr1 200 . A T . . . GT:PI:PG    0/1:3:0|1
chr1 300 . A T . . . GT:PI:PG    1/0:4:1|0
chr1 400 . A T . . . GT:PI:PG    0/1:4:1|0

where first two lines have same "PI" (i.e 3) stating they are phased with each other. Same with the later two lines which have another "PI" (i.e 4).


text:
Error in colSums(assay(sce)) : 'x' must be an array of at least two dimensions
I am having a SingleCellExperiment (sce) object:

sce <- SingleCellExperiment(list(counts=UMI_count))

which has two dimensions:
dim(sce)
[1] 46761 68027

When I am trying to get column sums for all rows the following does not work:
colSums(assay(sce), na.rm = FALSE, dims = 1)

Giving me an error:

Error in colSums(assay(sce)) : 'x' must be an array of at least two dimensions

Apparently the issue was already described:
https://support.bioconductor.org/p/100954/
However, I can not find the solution. Is there another way of getting sums of columns of either SingleCellExperiment object or just dgcMatrix?

Are you trying to get the sum of the rows or the columns?

Did you check what does assay(sce) return? Which is the output of dim(assay(sce)) ?

Ultimately I just implemented it manually with the precursor to SingleCellExperiment matrix (UMI_count which is dgcMatrix type) object in the following way:
# Removing genes with 0 - counts in all of the columns
num_removed = 0
for(i in 1:nrow(UMI_count)) {
  idx = i - num_removed
  if(sum(UMI_count[idx,]) < 1) {
    UMI_count <- UMI_count[-idx,]
    num_removed = num_removed + 1
    print(i)
  }
}

sce <- SingleCellExperiment(list(counts=UMI_count))

I am getting sums of the columns for each row and then filtering out rows with the sum < 1. So, sum function for the slices - rows (numeric type) is working.


text:
R: 'Matrix' can not be unloaded, but 'writeMM' method not found
I want to save a sparse and very large dgcMatrix onto disc. I read that it could be done with writeMM method. So, when I am trying to do that:
writeMM(UMI_count, "/home/gene_count_filtered/filtered_dataset")

Rstudio is giving me an error:
Error in writeMM(UMI_count, "/home/gene_count_filtered/filtered_dataset") :
could not find function "writeMM"

Ok, I read that this function is in the Matrix package. I installed it, but when I try loading it:
library(Matrix)

I am getting another error due to other libraries using Matrix already on the lower level, as I understand:
Error in value[[3L]](cond) :
Package ‘Matrix’ version 1.2.11 cannot be unloaded:
Error in unloadNamespace(package) : namespace ‘Matrix’ is imported by
‘scran’, ‘scater’, ‘SummarizedExperiment’ so cannot be unloaded

I also tried justwrite.matrix(UMI_count, "/home/nikita/filtered_dataset"), but got an error:

Cholmod error 'out of memory' at file ../Core/cholmod_memory.c, line 147

What could be the solution to the problem? I just need to somehow save a very large dgcMatrix on disc.

Which commands did you use to install Matrix? Did you use a fresh (with no environment)  session?

Nope, but I am having this object and I need to save it. I do not want to restart 10h calculation

you can install the package in a new instance of R

But how am I going to save the matrix of the current session? That is the issue

You open in terminal or from the menu another R session. From there you install the library, then go back to the terminal where you have the data, and try to save. If you find any problem, report the whole commands and output here.

Restarting R should fix this. There’s probably no way around this (you could try to manually unload all namespaces that use the Matrix package).

I was skeptical about it but restarting fixed everything.

Ultimately I just restarted Rstudio as Konrad Rudolph suggested, imported Matrix library, and then writeMM method worked just fine. However, rownames were converted into just numbers, but I needed to keep the gene names instead as it was in dgcMatrix, so after some search I found a different way of doing that from here:
http://www.sthda.com/english/wiki/saving-data-into-r-data-format-rds-and-rdata
So, I just did save() (and then, to load back - load()) and the rownames are maintained like I wanted it to be. So, storing it now as .RData file.


text:
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.

In my experience it's helpful to explore the data (ie. make some heatmaps etc) and look for trends that are associated with low molecule counts.  For example, cells with low umi's may look weird or cluster separately.  You can use these to guide your cutoff.

@GWW can you elaborate more on that? Do you mean for example to plot the two principal components and see where the different cells cluster? And maybe if I have a cluster composed by only cells with limited library size, exclude that cluster? Otherwise, which plots will you use? Thanks

I suggest a heatmap of the most variable genes.  PCA may be harder to see.

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).


text:
Installing RnBeads via bioconductor - .onLoad failed in loadNamespace(), call: NULL
I am trying to install RnBeads from bioconductor but the installation of its dependancy TxDb.Hsapiens.UCSC.hg19.knownGene fails. The error is:
biocLite("RnBeads")
BioC_mirror: https://bioconductor.org
Using Bioconductor 3.6 (BiocInstaller 1.28.0), R 3.4.3 (2017-11-30).
Installing package(s) ‘RnBeads’
also installing the dependencies ‘FDb.InfiniumMethylation.hg19’, ‘methylumi’

trying URL 'https://bioconductor.org/packages/3.6/data/annotation/src/contrib/FDb.InfiniumMethylation.hg19_2.2.0.tar.gz'
Content type 'application/x-gzip' length 38996322 bytes (37.2 MB)
==================================================
downloaded 37.2 MB

trying URL 'https://bioconductor.org/packages/3.6/bioc/src/contrib/methylumi_2.24.1.tar.gz'
Content type 'application/x-gzip' length 6178909 bytes (5.9 MB)
==================================================
downloaded 5.9 MB

trying URL 'https://bioconductor.org/packages/3.6/bioc/src/contrib/RnBeads_1.10.2.tar.gz'
Content type 'application/x-gzip' length 6437974 bytes (6.1 MB)
==================================================
downloaded 6.1 MB

* installing *source* package ‘FDb.InfiniumMethylation.hg19’ ...
** R
** data
** demo
** inst
** preparing package for lazy loading
Error: package or namespace load failed for ‘TxDb.Hsapiens.UCSC.hg19.knownGene’:
 .onLoad failed in loadNamespace() for 'TxDb.Hsapiens.UCSC.hg19.knownGene', details:
  call: NULL
  error: evaluation nested too deeply: infinite recursion / options(expressions=)?
Error : package ‘TxDb.Hsapiens.UCSC.hg19.knownGene’ could not be loaded
ERROR: lazy loading failed for package ‘FDb.InfiniumMethylation.hg19’

The downloaded source packages are in
    ‘/tmp/RtmpY8Yf4l/downloaded_packages’
installation path not writeable, unable to update packages: doMC, doSNOW,
  foreach, iterators, snow, boot, class, cluster, codetools, foreign,
  KernSmooth, lattice, MASS, Matrix, mgcv, nlme, nnet, rpart, spatial, survival

As suggested by the error message, I tried setting options(expressions=500000), which unfortunately did not change the error message. An important note may be that I am under network mount and don't have administrative access rights.
Edit
I don't have a .Rprofile file.
Note about installation path not writeable:
Typing .libPaths() reveals that my home R folder is at the first position (and it is writable). The other paths ("/usr/local/lib/R/site-library"                          , "/usr/lib/R/site-library" and "/usr/lib/R/library") are not writable without sudo.
I also tried to delete my home R folder and to reinstall it. The reinstallation of all packages worked for all packages except for the three listed above (same error message). I also tried biocLite("TxDb.Hsapiens.UCSC.hg19.knownGene", lib.loc=.libPaths()[1], lib=.libPaths()[1]) to specifically address my home R directory, but the same error is showing up (but without the installation path not writeable part).
Here is my sessionInfo():
R version 3.4.3 (2017-11-30)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.4 LTS

Matrix products: default
BLAS: /usr/lib/openblas-base/libopenblas.so.0
LAPACK: /usr/lib/lapack/liblapack.so.3.0

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

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

other attached packages:
[1] BiocInstaller_1.28.0

loaded via a namespace (and not attached):
[1] compiler_3.4.3 tools_3.4.3

Thank you very much!

It  might be better then to raise the error to support.bioconductor.org to let the maintainers know about it.

Ok, I will do. If the error is fixable from user side, I will report the answer here.

I have reported the error on support.bioconductor.org but we were not able to solve the problem. It seems that some dependencies are making problems, but as reinstalling my home R folder did not resolve the issue, I have come up with the following workaround:
I downloaded the latest R version from my favorite CRAN mirror, opened an R session with pathToNewR/bin/R from my terminal and installed RnBeads there. This works fine.
Note: I had to make sure that pathToNewR/library is the only entry in .libPaths(), as per default, my home R folder was also contained in .libPaths().
As a consequence, I am independent from the (maybe broken) R version and packages on my server and can finally work with RnBeads :)

If I come up with another solution that works fine with my server R version, I will post an update here.


text:
How to compare transcriptomic profiles of two cell types (single cell RNA-seq)?
I found this interesting Single RNA-seq data set in GEO, but I am not sure how to analyze it appropriately.
They have deposited transcriptomic profiles of human and mouse pancreatic islets (pancreatic cells: Beta cells, Delta, etc).  The problem I see is that the different pancreatic cell types are not in equal numbers: 872 Beta cells and 214 Delta cells.
I want to see whether the expression of a set of genes are correlated between the two cell types. I don't think that I can visualize the correlation between the expression of the Beta and Delta cells on a scatter plot if they are unequal in number. Is there a workaround to include all the cells in the analysis, without losing good quality transcriptome of beta/delta cell?
What I would like to do is to compare the expression of two pancreatic cell types (Beta v.s. Delta cells) using a scatter plot. I have a set of genes that I am interested in, and would like to see their expression profile for beta and delta cells. As far as I know, it is only possible to visualize this using a scatterplot (X:Beta vs Y:Delta) where the number of isolated beta cells is the same as the number of delta cells.
Question: Given unequal number of isolated panacriatic cells, what would be an appropriate way to compare the expression profiles? Should I just ignore the extra ones?
Any idea how to approch this problem?
Thanks,

What is wrong with a traditional differential expression analysis? How big is the difference in number between Beta and Delta cells?

@Llopis Quite large ... the difference between beta and delta is 658 cells.

Do you have a particular set of genes in mind or are you hoping to go fishing for one? It's also unclear why you can't do visualizations with unequal sample sizes.

Any clarifications (such as what you have just written) should be put in the question, rather than the comment.

What makes you think you wouldn't be able to visualise the correlation for unequal numbers?

@MEhsan You can randomly pick the same amount of delta cells than of beta cells and bootstrap to find the real correlation (or plot the scatter plot)

@Llopis I'm not familiar of how to use bootstrapping and randomly sample the cells.  Could you give me an example?  Is there a way to do this in Python ...?

It is "just" a for loop where you pick different samples, I have never done it in Python so you'll need to search in the web.

You are looking for a way to compare expression profiles between cell types, and your data is counts of genes expressed per cell. Your issue is that the data is highly dimensional. It has as many dimensions as you have genes. To be able find useful information, you must reduce these dimensions using a statistical technique such as Principal COmponents Analysis (PCA), t-Distributed Stochastic Neighbor Embedding (t-SNE), or Multidimensional Scaling (MDS). The result of these analyses is a reduction of your thousands of dimensions to vectors (eigenvectors in PCA) of variance which you can plot against each other in a scatter plot.
I have been working with Single Cell RNA-seq data and have been exploring the difference in gene expression in a certain type of cells between timepoints. Seurat (R package) is a very powerful tool for scRNA-Seq data analysis.
Here is a a tSNE of my data. I have added the timepoint labels, but you can see that after performing the dimensionality reduction the two different timepoints cluster by their gene expression profiles. I can then use Seurat to see which specific genes are differentially expressed between the clusters and a whole bunch of other fun things.

I just add the module to run t-SNE in python, since OP seems to want to use a python module: http://scikit-learn.org/0.15/modules/generated/sklearn.manifold.TSNE.html

Scatterplots are used to represent multiple dimensions of the same piece of data. A scatterplot would be used in this case to display all genes in your set on a single plot, but I would have some difficulty interpreting such a plot. If you want to look separately at each gene, then a density curve (or histogram) would be more appropriate.
Your problem is similar to common problems in flow cytometry, and I don't think there's any reason why similar approaches can't be used. Here's a histogram / density plot showing two different types of cells (negative control, and treatment condition), indicating the difference in fluorescence of a particular marker (somewhat analogous to gene expression) across all cells of each type. Note that the number of cells are different between the negative control and treatment:

This image was taken from here, which also demonstrates a 2D scatter plot (where each point is a cell, rather than a gene).

seurat3 may help. You can compare the marker genes from the 2 cell types or treat one cell type as a control. Please see the example analysis from: https://satijalab.org/seurat/v3.0/immune_alignment.html


text:
What software/tool can I use to visualise a multiple sequence alignment in the xmfa format as produced by command-line tool Mauve?
What software/tool can I use to visualise a multiple sequence alignment in the xmfa format as produced by command-line tool Mauve?
I am particularly interested in seeing the single nucleotide polymorphisms on the screen.
I am doing a multiple sequence alignment of about a dozen whole genomes, all in a single contig.
Link to Mauve

What sort of alignment? How may sequences? Are you thinking reads against genome, or is this multiple full length genomic sequences aligned against each other? Since you mention mauve, should we assume whole genome alignments?

Whole genome alignments of about a dozen full bacterial chromosomes.

Typically I use the Mauve viewer to visualize the alignments.  Are there things you find lacking in that viewer?
There are multiple threads about converting xmfa to other formats

I would like to see the sequences more zoomed in, at the base level. Do you know if it's possible with the Mauve viewer?

I don't think so.  Mauve is a designed to find large segments of homology and is not really for per base alignments

I supervised a student a while ago, and she made a tool to focus on (i.e., visualize) SNPs in multiple alignments. The tool is called ADOMA, and is quite simple. It uses clustal for alignment, so I am not sure it if it may work for your Mauve data. If your data are xmfa, you'll probably have to make fasta of it first.
So I am not sure if this is a solution, but maybe the tool can be of use for your research.

Can ADOMA be used for whole genome alignments, e.g. bacterial genomes? Or is just for gene-size alignments?

Sorry I am not sure. The tool was made by a very good (independent working) student. It uses clustal, and if I recall correctly, clustal can also take aligned data already. So I guess your aligned data can maybe be formatted to the right input for clustal? And then adoma could in theory show only SNPs in the alignment. I think whole genome data would make very large output files, so I am not sure if it would be feasible.


text:
Map domain names from UniProt bed files to domain accessions
I want to get a bed file mapping human protein domains to the human genome. UniProt actually offers such a thing here. The problem, however, is that the file doesn't include any kind of domain accession, so I have no way of knowing exactly what domain this is. For example:
$ grep P30443 UP000005640_9606_domain.bed | sed -E 's/\t+/ /g'
chr6 29933668 29933929 P30443 0 + 29933668 29933929 153,153,255 1 261 0 . P209-T295; Ig-like C1-type
chr6 29936281 29936542 P30443 0 + 29936281 29936542 153,153,255 1 261 0 . P209-T295; Ig-like C1-type

This tells me there's a "Ig-like C1-type" domain at positions 209-295 of the protein. But what is this domain?
The UniProt page for this protein doesn't add any more information:

I thought they might be using the official description of the domain from PFam. However, there is no such domain in PFam as far as I can tell. I downloaded the full PFam annotation of SwissProt from PFam's FTP site and searched for P30443:
$ zgrep -A4 P30443 swisspfam.gz | sed -E 's/  */ /g; s/[-=_]+//g'
>1A01HUMAN || P30443.1 365 a.a.
C1set 1  (5782) PF07654.14 Immunoglobulin C1-set domain 210290
MHCI 1  (1847) PF00129.17 Class I Histocompatibility antigen, domains alpha 1 and 2 25203
MHCIC 1  (439) PF06623.10 MHCI Cterminus 337364

This tells me there is a "Immunoglobulin C1-set domain" at 210-290. So both the coordinates and the domain "name" are different.
I also checked the Pfam-A.regions.uniprot.tsv.gz file available here. This one had:
$ zgrep P30443 Pfam-A.regions.uniprot.tsv.gz | sed -E 's/\t+/ /g'
P30443 1 8667AFF3F06C4932 d5f7ef1d93bd6e2e0ca1f7b3477ffb67 PF07654 210 290

This gives me a PFam accession (PF07654) but, once more, I can't use this to go back to the bed file because the coordinates are different.
I hoped I could figure this out using the SwissProt flat files, but no luck there either. These are the relevant lines from the entry for P30443:
ID   1A01_HUMAN              Reviewed;         365 AA.
AC   P30443; O77964; O78171; Q9MYA3; Q9TP25; Q9TQP5;
[. . .]
DR   InterPro; IPR021157; Cyt_c1_TM_anchor_C.
DR   InterPro; IPR007110; Ig-like_dom.
DR   InterPro; IPR036179; Ig-like_dom_sf.
DR   InterPro; IPR013783; Ig-like_fold.
DR   InterPro; IPR003006; Ig/MHC_CS.
DR   InterPro; IPR003597; Ig_C1-set.
DR   InterPro; IPR011161; MHC_I-like_Ag-recog.
DR   InterPro; IPR037055; MHC_I-like_Ag-recog_sf.
DR   InterPro; IPR011162; MHC_I/II-like_Ag-recog.
DR   InterPro; IPR001039; MHC_I_a_a1/a2.
DR   InterPro; IPR010579; MHC_I_a_C.
DR   Pfam; PF07654; C1-set; 1.
DR   Pfam; PF00129; MHC_I; 1.
DR   Pfam; PF06623; MHC_I_C; 1.
DR   PRINTS; PR01638; MHCCLASSI.
DR   SMART; SM00407; IGc1; 1.
DR   SUPFAM; SSF48726; SSF48726; 1.
DR   SUPFAM; SSF54452; SSF54452; 1.
DR   PROSITE; PS50835; IG_LIKE; 1.
DR   PROSITE; PS00290; IG_MHC; 1.
FT   DOMAIN      209    295       Ig-like C1-type.

So, while it does have this domain, uses the same "name" and coordinates I had in the bed file ("Ig-like C1-type"), and it has the relevant PFam accession (PF07654), there's no way to link the FT line with "Ig-like C1-type" to the relevant PFam accession.
I even tried using ProSite and downloaded prosite.dat from here. This one has yet another "name" for the domain:
AC   PS50835;
DT   01-MAY-2003 CREATED; 10-MAY-2017 DATA UPDATE; 20-DEC-2017 INFO UPDATE.
DE   Ig-like domain profile.

Finally, I had hoped I could use the coordinates given in the bed file to map but since those aren't the same as those I find in the various PFam files, I'm stumped.
Is there any way to map the domain names given in the UniProt bed files to specific domains (either PFam, or PROSITE or anything else) so that I can know exactly what domain is annotated there?
Update: After jumping through various hoops using the files mentioned above, I have managed to map only 701 domain names which leaves me with 2764 unmapped.


text:
Using column 2 of one file to match with two columns of another file, and append
I have file 1 like following:

1   15776220    15776240    GTGACCAGCAGGTGTCTCTG    16855676    16855696    CTGTCCAGCAGAGGGCGGTG

And file 2 as following

1   15776231    2   5008    G:5002  A:6 1   16855677    2   5008    A:5003  C:5

I am trying to write a code(but failing!) to do the following:
if $2 of file2 comes in between $2 and $3 of file1 (between the interval of $2 and $3) and if $8 of file2 comes in between $5 and $6 of file1 (between the interval of $5 and $6), then the output would be $1,$2,$3,$4 from file1 and $2,$3,$4,$5,$6 from file2 and $5,$6,$7 from file1 and $8,$9,$10,$11,$12 from file 2, all in one line. So something like this:
1   15776220    15776240    GTGACCAGCAGGTGTCTCTG    15776231    2   5008    G:5002  A:6 16855676    16855696    CTGTCCAGCAGAGGGCGGTG    16855677    2   5008    A:5003  C:5

Something like bedtools does but I am not able to use bedtools one this since I don't have two columns in file2. Is there any way to make bedtools work on this?
(There could be many lines from file2 intersecting between the same line in file1,so that line could be repeated)

I see you took my advice, welcome to [bioinformatics.se]! How big are the files? Are they small enough to fit in RAM? Can we load one of the files in memory and then process the other?

Yes they can be worked on in the RAM itself, big but not huge.

Pipe a modified form of the second file to BEDOPS bedmap and the first file, then pipe that result to cut out desired columns:
$ awk -v OFS='\t' '{ print $1, $2, ($2+1), $3, $4, $5, $6, $7, $8, $9, $10; }' second.txt | bedmap --echo --echo-map first.bed - | cut -f1-4,8,9,11- - > answer.bed

This should run pretty quickly and use very little memory, as it uses Unix streams on (presumably) sorted BED data.
You'll want to test the column indices in my answer by running this on your test inputs, but that's an easy adjustment for the larger dataset, in case I missed a column.
If you're not sure about sort order, pipe to sort-bed:
$ sort-bed first.unsorted.bed > first.bed
$ awk -v OFS='\t' '{ print $1, $2, ($2+1), $3, $4, $5, $6, $7, $8, $9, $10; }' second.txt | sort-bed - | bedmap --echo --echo-map first.bed - | cut -f1-4,8,9,11- - > answer.bed

Use streams where you can. It's a different way of thinking but can pay huge dividends.

Nice! What is the ($2+1) doing in your last awk command? I mean, I know what it does, it prints the value of the second field plus one, but why?

This creates a half-open interval that can be used for set operations. I'm assuming that both files are on the same coordinate system. If not, that can be adjusted with awk.

Hey thanks for the script, but it doesn't give me any intersection. It only repeats $2 $3 and $5 $6  with the sequence in some lines, and I get no values from my second file in my answer

Here's one approach in Perl. Note that I have not tested this extensively, and I am not confident it works. Please test first and let me know if there are problems.
#!/usr/bin/env perl
use strict;
use warnings;
use 5.10.0;
my (@file1,%file2);
open(my $f1, '-|',"sort -nk2 $ARGV[0]");
while (<$f1>) {
  chomp;
  my @fields = split(/\s+/);
  push @file1, [ @fields ];
}
close($f1);

open(my $f2, "sort -nk2 $ARGV[1] |");
line:while (<$f2>) {
  chomp;
  my @fields = split(/\s+/);
  foreach my $file1Line (@file1) {
    my $file1Start1 = $file1Line->[1];
    my $file1End1 = $file1Line->[2];
    my $file1Start2 = $file1Line->[4];
    my $file1End2 = $file1Line->[5];
    my $file2Start1 = $fields[1];
    my $file2End1 = $fields[2];
    my $file2Start2 = $fields[7];
    my $file2End2 = $fields[8];
    my @file1Fields = @$file1Line;
    if ($file1Start1 < $file2Start1 &&
        $file1End1 > $file2End1 &&
        $file1Start2 < $file2Start2 &&
        $file1End2 > $file2End2) {
      say join "\t", @{$file1Line}[0..3], @fields[1..5], @{$file1Line}[4..6], @fields[7..11];
    }
    elsif ($file1Start1 > $file2Start1) {
      last;
    }
  }
}
close($f2);

Save the script as foo.pl somewhere in your $PATH, make it executable and run as:
$ foo.pl file1.bed file2.bed
1   15776220    15776240    GTGACCAGCAGGTGTCTCTG    15776231    2   5008    G:5002  A:6 16855676    16855696    CTGTCCAGCAGAGGGCGGTG    16855677    2   5008    A:5003  C:5

Hey thank you for the script. I think it almost works but there are problems. Shell gives a few Argument errors ("Argument "end" isn't numeric in numeric gt") for a lot of lines and "Use of uninitialized value $file2Start2 in numeric lt" for a lot of lines as well.
However it does output a file, but with incorrect intersections and the number of lines in the file is 30 times more than the first file, which shouldn't be.
The files are only 90kb and 30kb in size. Would you mind taking a closer look if I e-mailed them to you or something?

@MaharshiChakravortee this isn't a shell script. It's a perl script. Run it as I say in my answer, but making it executable and placing it in your PATH, or run it like this: perl foo.pl file1.bed file2.bed . If this still fails, ping me (@terdon) in chat.


text:
How to remove all BAM read groups from all reads (not just the header)?
I have problem with one my BAMs---it appears to have invalid read groups.
Normally when I have such a problem, I remove all the read groups from the BAM header as follows:
samtools view -H your.bam | grep -v "^@RG" | samtools reheader - original.bam > new.bam

and then I would execute AddOrReplaceReadGroups from Picard:
http://broadinstitute.github.io/picard/command-line-overview.html#AddOrReplaceReadGroups
I have tried this, and there must still be reads with a read group in the optional tags RG:Z:.
How to do I definitively remove all read groups from a BAM, including all of the reads?

You might have to manually strip those auxiliary tags off:
samtools view -h your.bam | grep -v "^@RG" | sed "s/\tRG:Z:[^\t]*//" | samtools view -bo your_fixed.bam -

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


text:
replacing SNPs with missing calls with a specific string
I have a big file containing 27 columns and nearly 6 million rows. The following is a little example of my file
head data
0.65   0.722222   1.0      0.75     0
0.35   0.277778   0.0      0.25     0
0      0.666667   0.75     0.5      0.5625
0      0.333333   0.25     0.5      0.4375

Rows are SNPs and columns are individuals. and I have "2 rows per each SNP" (one for reference allele and the other for alternative allele). In the example above I showed data for 2 SNPs (rows 1 and 2 correspond to SNP1 and rows 3 and 4 correspond to SNP2). I want to check if for each SNP, both ref and alt allele frequencies are 0, replace them with 9. this is my desired output:
head desired
0.65   0.722222   1.0      0.75     9
0.35   0.277778   0.0      0.25     9
9      0.666667   0.75     0.5      0.5625
9      0.333333   0.25     0.5      0.4375

Do you only want to replace a field with 9 if the same field of the other line is 0 or if any of them are? For example, if you have a line with 0.65   0.722222   1.0      0.75     0 and the corresponding line is 0   0.722222   1.0      0.75     0.3 (so field 1 is 0 in one line and the last field is 0 in the other), what should happen? Should both 0s be made 9s? Neither?

A possible approach in plain python:
#!/usr/bin/env python3

import sys

def both_zero(freq1, freq2):
    return freq1 == "0" and freq2 == "0"

with open(sys.argv[1], "r") as snp_file:
    for line1 in snp_file:
        # Consume one extra line to get the frequencies
        # for the alternative allele
        line2 = snp_file.readline()
        # "zip" the frequencies in pairs
        freq_pairs = zip(line1.strip().split(), line2.strip().split())
        recoded_pairs = [
            ("9", "9") if both_zero(freq1, freq2) else (freq1, freq2)
            for (freq1, freq2) in freq_pairs]
        # "unzip" the recoded frequencies into separate objects
        # (see https://stackoverflow.com/a/19343/1878788)
        (recoded1, recoded2) = zip(*recoded_pairs)
        print(*recoded1, sep="\t")
        print(*recoded2, sep="\t")

Calling it with your example snp frequencies file as command-line argument:
$ ./recode_absent.py snps.txt
0.65    0.722222        1.0     0.75    9
0.35    0.277778        0.0     0.25    9
9       0.666667        0.75    0.5     0.5625
9       0.333333        0.25    0.5     0.4375

thanks. I just tried your python code, but I got a syntax error:  File "./replace.py", line 21
print(*recoded1, sep="\t")
^
SyntaxError: invalid syntax

@Anna1364 Maybe you are using python2 and not python3. Either try to add from __future__ import print_function at the top of the script, or try to execute the script with python3


text:
I am getting this error when running a test data set for Lyve-SET on RHEL 7. Is this an MPileup error or something else?
This is the error that I am getting:
launch_varscan.pl: Creating a pileup /home/user/lyve-SET-1.1.4f/lambda/tmp/NC001416.fasta.wgsim.fastq.gz-reference.sorted.bam.mpileup
removed ‘/home/user/lyve-SET-1.1.4f/lambda/tmp/NC001416.fasta.wgsim.fastq.gz-reference.sorted.bam.1239.mpileup’
launch_varscan.pl: Running mpileup:
  echo "NC_001416.1
" | xargs -P 1 -n 1 -I {} sh -c 'echo "MPileup on {}" >&2; samtools mpileup -f lambda/reference/reference.fasta
-Q 0 -B  --region "{}" /home/user/lyve-SET-1.1.4f/lambda/bam/NC001416.fasta.wgsim.fastq.gz-reference.sorted.bam >
/home/user/lyve-SET-1.1.4f/lambda/tmp/NC001416.fasta.wgsim.fastq.gz-reference.sorted.bam.$$.mpileup'
MPileup on NC_001416.1
mpileup: unrecognized option '--region'
open: No such file or directory
[mpileup] failed to open NC_001416.1: No such file or directory
launch_varscan.pl: main::mpileup: ERROR with xargs and samtools mpileup
Stopped at /home/arahat/lyve-SET-1.1.4f/scripts/launch_varscan.pl line 106.
QSUB ERROR
31488
launch_set.pl: Schedule::SGELK::command: ERROR with command: Inappropriate ioctl for device
  /usr/bin/perl /home/user/lyve-SET-1.1.4f/lambda/log/SGELK/qsub.2705.pl
Stopped at /home/user/lyve-SET-1.1.4f/scripts/../lib/Schedule/SGELK.pm line 661.
set_test.pl: ERROR: died while trying to run COMMAND
  launch_set.pl lambda --numcpus 1 --numnodes 20
set_test.pl: See if you have the correct software installed for this command.
6400 at /home/arahat/lyve-SET-1.1.4f/scripts/set_test.pl line 94.
make: *** [test] Error 25

I think the main problem is the "mpileup: unrecognized option --region", which seems to be blank. In the beginning of the block, it says creating pileup and then it immediately removes the file, probably to create the pipeline file.
I cannot change the code for Lyve-SET pipelining or create shell scripts for the pipeline because I am only configuring my server for use with the genome sequencing pipeline. Would this be a problem with the data set or some tool that I need to install? Any insight would be appreciated.

What version of samtools mpileup to you have installed?  What version of varscan are you using?  Did you install the two tools separately or did varscan come with samtools?  I have version 1.3.1 and there is a --region option for mpileup.

Have you tried running samtools mpileup -f lambda/reference/reference.fasta  -Q 0 -B  --region NC_001416.1 /home/user/lyve-SET-1.1.4f/lambda/bam/NC001416.fasta.wgsim.fastq.gz-reference.sorted.bam  directly to see if the command works outside the context of varscan?

I installed samtools version 0.1.19 using yum install. Varscan was installed separately through make install from Lyve-SET. I did not try to run samtools mpileup separately because all the commands are implemented through Lyve-SET, which is how it will be run by the user. I will try reinstalling varscan and running the samtools mpileup command separately. Thank you for the help.

I would suggest installing a more recent version of samtools.  That version is from 2013 is likely is not compatible with the version of varscan you installed

Yup. I noticed that yum patches but does not update the packages so I am manually installing the newest versions of all the tools specified in Lyve-SET. Thank you for bringing this to my attention. Hopefully it works after I reinstall.

The samtools v0.1.19 manual confirms that there is no --region option.  In that version there is only the -r argument to pass in the targeted region.  If you update to a newer version of samtools this should execute correctly

It worked after I reinstalled all the tools. The main problem was the mpileup error where --region was not an option. Cannot trust some of these yum packages as they are much older versions. Thank you so much for helping me with this issue.


text:
Scaling by linear regression against the number of reads
I am trying to build the preprocessing pipeline presented in The Tabula Muris Consortium et al. (pp).
It is a pipeline to preprocess single-cell sequencing data. There is one step that is not clear:

Counts were log-normalized (log(1 + counts per N)), then scaled by linear regression against the number of reads (or UMIs), the percent of reads mapping to Rn45s, and the percent of reads to ribosomal genes.

I understand the first part (I assume that log in this context is log2), but I need help on understanding how to scale by linear regression against the number of reads, the percent of reads mapping to Rn45s, and the percent of reads to ribosomal genes.

Have you contacted the paper authors? They would be able to give you a better answer, and your comments/questions would help them improve the paper when it is properly published.

@gringer You are right, I posted here because I thought there was some standard procedure.

The linear regression is straight-forward enough, but regressing against only three values (if I'm interpreting that correctly) is a recipe for disaster (it sounds vastly worse than even ERCC spike-ins...and those aren't exactly ideal).

@DevonRyan ignoring by now the small number of regression values, how do you scale using the linear regression? I am unable to see how the linear regression is used in this case. Can you elaborate more (maybe with an answer)?

Did you figure how did they preprocess the tabula muris data?

@yuqi_yuqi see https://bioinformatics.stackexchange.com/a/3225/1771

I don't know if this question has been solved already, but what they try to do is equalize the depth of sequencing for each cell. Therefore, they scale for the total number of reads. If you regress out (via linear or negative binomial regression) the differences in the number of reads per cell, you end up with cells that have been sequenced with the same depth.
In my opinion they use the ribosomal genes in the same way. They are considered to be a kind of housekeeping genes, which you can use to equalize the sequencing depth.
And yes, the log in this context is log2, used to get fold-change values instead of counts.


text:
What are doublets in single cell RNA-seq data?
I am reading The Tabula Muris Consortium et al. (pp).

In some organs, cells with more than 2 million reads were also excluded as a conservative measure to avoid doublets.

How exactly is a “doublet” defined? For example, is doublet a set of cells sequenced as a single cell (and so, the number of transcripts is double)?

Care to explain downvotes?

I haven't downvoted you but this question (to my taste) doesn't show any research effort, I think this kind of downvotes can be prevented by showing some effort to solve yourself the question, for example, showing which articles have you read about doublets (if they exists) or what kind of search have you done to find more information, even if it is unsuccessful.

Is doublet a set of cells sequenced as a single cell?

Yes. Depending on the method of single cell sequencing it may be more or less likely for groups of cells to be captured and barcoded with the same "unique" barcode. This is more likely in split-pool RNA sequencing (e.g. SPLiT-seq), and less likely in cell-capture RNA sequencing (e.g. Fluidigm C1). Doublets can also be created through physical / experimental processes (e.g. from tissues that were not completely dissociated).
Bear in mind that detecting doublets is not as simple as counting for doubling of transcript expression, because the expression profiles are different in different cells. That's why single-cell sequencing is useful in the first place.

“less likely in cell-capture RNA sequencing” — Except that the C1 used to have a notorious flaw that messed up its doublet detection when they were z-stacked.

"Doublet" is commonly used to describe a droplet in droplet-based sequencing that has captured atleast 2 cells. 10x states their doublet rate to be 0.8% per 1000 cells:

There is a tradeoff between targeted cell capture and doublet rate and lab protocols are usually optimized for a certain expected doublet rate.
One should expect to see close to doubling of transcripts but it is not trivial to infer this computationally.


text:
Can the canu assembler output a fastq file of the final assembly just like HGAP4?
I have assembled some genome from Sequel PacBio data both with HGAP4 on the SMRT Link interface and using canu on the command line. The HGAP4 assembler outputs a fastq file of the final assembly such that I can see how reliable the sequence is along the genome. On the other hand, canu does no such thing, which means I do not know how much I can trust each nucleotide in the assembly? Is there a way canu can output the fastq file as well or at least some depth information at each base?

I've created an issue on the Canu github repository for this. I'm not aware of any existing functionality to output FASTQ files, but think that this would be a useful feature to have.
It would be possible to create something like this by aligning the original reads to the genome, generating a pileup from that alignment, then determining the sequence likelihood based on that pileup. Unfortunately, I don't have any scripts that I've written to do that... yet.


text:
Biopython - extracting query coverage from XML Blast output
I need to extract (or to calculate precisely) the query coverage from Blast output. I am using NCBIWWW.qblast to query Blast in XML format.
Is it possible to know the query coverage parsing the XML output? If yes, how?
result_handle_xml = NCBIWWW.qblast("blastp","pdb",my_query.seq,format_type='XML')

here it suggest to use '-6 qcovs', but only for BLAST command line version and only for tabular format
This is the fasta file that I open:
>sp|P02489|CRYAA_HUMAN Alpha-crystallin A chain OS=Homo sapiens
GN=CRYAA PE=1 SV=2
MDVTIQHPWFKRTLGPFYPSRLFDQFFGEGLFEYDLLPFLSSTISPYYRQSLFRTVLDSG
ISEVRSDRDKFVIFLDVKHFSPEDLTVKVQDDFVEIHGKHNERQDDHGYISREFHRRYRL
PSNVDQSALSCSLSADGMLTFCGPKIQTGLDATHAERAIPVSREEKPTSAPSS

Please [edit] your question and include an example of the file you are trying to parse. We can't really help you parse data you don't show.

I see you edited, but you still haven't shown us the file you are trying to parse. Is the information even there? Please add an example so we know what you're working with.

According to this, the answer is 'no'.

@AndreaSpinelli please consider the answers provided and whether these answered your question. If so please remember to "accept" an answer, it helps with site stats

I had the same problem. I'm pretty sure you can calculate it yourself using the following:

The length of your sequence (seq_len)

The value from the <Hsp_query-from> tag from the hit (<Hsp_query-from>)

The value from the <Hsp_query-to> tag from the hit (<Hsp_query-to>)

It can be calculated by:
query_cover = (<Hsp_query-to> - <Hsp_query-from> + 1) / seq_len * 100

Please note that I did not get this from any documentation. I figured it out by staring at the xml output returned from the Biopython NCBIWWW.qblast() function. This solution seems to match with the results I get by blasting the same sequence on the interactive web portal.

record = NCBIXML.read(result_handle_xml)
query_len = record.query_length
for alignment in record.alignments:
    _ = alignment.hsps.__next__()
    query_cover = (_.query_end - _.query_start + 1) / query_len * 100


text:
Inspection of gene expression in scRNA-seq data
I am running the data preprocessing pipeline for scRNA-seq data presented here.

3.8.6.1 Gene expression
In addition to removing cells with poor quality, it is usually a good idea to exclude genes where we suspect that technical artefacts may have skewed the results. Moreover, inspection of the gene expression profiles may provide insights about how the experimental procedures could be improved.
It is often instructive to consider the number of reads consumed by the top 50 expressed genes.

By analyzing my data I have the following most expressed 50 genes (evaluating the total number of transcripts across cells):
Original

Log-transformed

Since I'll be working on log-transformed data, I assumed that the distributions of the top 50 genes were not excessively skewed. However, as you can see in the first plot, there are some samples (cells) with very high number of transcripts for the top gene. Is it something I should worry about? What other insights can I get from these two figures?
Edit
Original analysis (figure attached) was:

The distributions are relatively flat indicating (but not guaranteeing!) good coverage of the full transcriptome of these cells. However, there are several spike-ins in the top 15 genes which suggests a greater dilution of the spike-ins may be preferrable if the experiment is to be repeated.

In an ideal world, ribosomal RNA (as seen in your top hit) should be excluded from samples prior to sequencing. Where this is not possible (i.e. in the data that have been presented to you), it would be a good idea to exclude ribosomal genes (and any other common contaminants) prior to doing further analysis (including normalisation). I believe this is the point the course instructors are trying to get across with that message box.
Eyeballing the results, you should also exclude any ERCC genes as well, as they are synthetic constructs and not from the target samples (although record their counts, as it may be informative for QC/normalisation purposes); they'll have odd distributions because of this.
ERCC sequences are typically added based on the concentration of DNA in a sample, rather than the number of cells. However with single-cell data (i.e. with known cell counts), it'd be possible to do a proper ERCC-based normalisation, assuming an equal amount were added for each cell. I suspect that knowing a bit more about the ERCC addition protocol would help in understanding why their distribution is so different.

Thanks! I noticed that the original comment for the analysis was not included in the question: I updated the question, I didn't get what they meant with flat distribution.

I forgot to ask, why remove ribosomal genes? I tried to find an explanation of this contamination, but with no results.

As you've noticed from your own analysis, the ribosomal genes have quite variable expression across cells. They're expressed everywhere, and quite difficult to completely remove from a sample (even with polyA selection), with success of removal depending on things like the amount of mechanical damage that the RNA has been exposed to. This high, variable expression throws out normalisation calculations, so it's better to remove them first.