RNAseq: Z score, Intensity, and Resources I'm very new to bioinformatics in general, and I'm trying to understand some basic concepts.

I have RNAseq data, and bioinformatics people tell me that intensities cannot be compared across patients. So there are all of these pipelines to compare intensities to Z scores--are those as simple as just plugging data into a bioconductor package?

It would be great to get some overview description of why Z values are important/why you can't compare intensities across patients, and/or some pointers towards any resources that I can read to learn more.

Thanks much for your time!

It depends on what test or analysis you want to do, whether you need intensities (expression values) or z-scores.

If you want to do statistical analysis, such as finding differentially expressed genes between groups of patients (e.g., with limma), you don't want to use z-scores. But you use normalized intensities (for microarrays) or start with raw counts for RNAseq. Then follow the user guide (limma has a great user guide, also for starters).

For visualization in heatmaps or for other clustering (e.g., k-means, fuzzy) it is useful to use z-scores. Z-scores are a form of transformation (scaling), where every genes is sort of "reset" to the mean of all samples, using also the standard deviation. If you want to know exactly what a z-score is, a simple google search can tell you the details.

In R you can use the scale function for z-score transformation. Be aware that the function works on columns though. Which means you have to transpose your matrix first if genes are in the rows, and then transpose them back after transformation.


How to identify the reference genome assembly used to create a .hic file? I have a bunch of .hic files of different species. But, I don't know which reference genome assembly (e.g. GRCh37 or GRCh38 in the case of human genome) the genomic co-ordinates of the .hic file correspond to. So I wonder if there is any way to figure out the version of the reference genome assembly from a .hic file.

Therefore, my question is: Does .hic file contain metadata such as the version of reference genome assembly? If yes, then how could I access it?

(I could obtain the information I need by digging into the methods sections of each paper corresponding to each .hic file. But it would be much better if I could get that right-away from the .hic file.).

Turns out HiC file does store such metadata.

straw tool from Aiden Lab contains a script (read_hic_header.py) that prints out the metadata from a given HiC file. Usage: read_hic_header.py <hic file or URL> [verbose]

In order to obtain the metadata in the form of a dictionary, one may consider using a function of the same tool that is currently available in my fork.

from straw import straw_module
straw_module.read_metadata(infile,verbose=False)

Args
    infile: str, path to the HiC file
    verbose: bool, print the metadata if True else return the metadata
Returns
    metadata: dict, the metadata is returned if verbose is False.
                Keys of the metadata:
                HiC version,
                Master index,
                Genome ID (str),
                Attribute dictionary (dict),
                Chromosomes (dict),
                chromosomes2length (dict),
                Base pair-delimited resolutions (list),
                Fragment-delimited resolutions (list).


featureCounts segmentation fault on Arch Linux I am encountering a segmentation fault when attempting to run featureCounts from subread-1.6.3 on even small test data. I installed featureCounts from SourceForge using the precompiled Linux binary.

Most existing help threads on this topic (1, 2, 3) show some splash screen output before the error, but I get a segmentation fault before any error is shown.

$ featureCounts -T 1 -p -t exon -g gene_id -a /data/sai/annots/musmus_cdna_all_ensembl92.gtf -o test.txt test.bam test2.bam
Segmentation fault (core dumped)

I resolved this problem by installing from source; I assume the precompiled binary was not compatible with the Arch Linux kernel.


Create identicon of a DNA sequence How can I create an identicon (a visual representation) of a DNA sequence?

Sequenticon is a Python library for generating identicons for DNA sequences. For instance, the sequence ATGGTGCA gets converted into the following icon:

A web interface is also provided at EGF CUBA (Collection of Useful Biological Apps): Render Sequenticons
Disclaimer: I'm the current maintainer of Sequenticon


Merge hundreds of small BAM files into a single BAM file I am working with over a million (long) reads, and aligning them to a large genome. I am considering running my alignment jobs in parallel, distributing horizontally across hundreds of nodes rather than trying to run a single job with dozens of cores.

I would like to merge the sorted BAM files together for further downstream analysis. What is the most efficient way to do so while maintaining a valid file header and taking advantage of the fact that the input bam files are already sorted?

samtools merge merged.bam *.bam is efficient enough since the input files are sorted. You can get a bit faster with sambamba and/or biobambam, but they're not typically already installed and IO quickly becomes a bottleneck anyway.


Get the mapping statistics of a single read from a BAM file I have a BAM file, and I have a read ID. What is the simplest way to get mapping statistics of that read in human-readable format?

E.g. I might want: % identity of aligned bases; number of insertions and deletions; number of bases aligned to the reference; etc.

Extra points for one-liners that give detailed summary statistics, and extra-extra points if you can extend the answer to output summaries for all reads longer than some predetermined length.

samtools stats seems to be able to do most of this, excluding the CIGAR-string parsing stuff (i.e. INDELs):

$ samtools view -h mapped.bam | grep -e '^@' -e 'readName' |
    samtools stats | grep '^SN' | cut -f 2-

raw total sequences:    2
filtered sequences: 0
sequences:  2
is sorted:  1
1st fragments:  2
last fragments: 0
reads mapped:   2
reads mapped and paired:    0   # paired-end technology bit set + both mates mapped
reads unmapped: 0
reads properly paired:  0   # proper-pair bit set
reads paired:   0   # paired-end technology bit set
reads duplicated:   0   # PCR or optical duplicate bit set
reads MQ0:  0   # mapped and MQ=0
reads QC failed:    0
non-primary alignments: 0
total length:   13632   # ignores clipping
bases mapped:   13632   # ignores clipping
bases mapped (cigar):   13502   # more accurate
bases trimmed:  0
bases duplicated:   0
mismatches: 3670    # from NM fields
error rate: 2.718116e-01    # mismatches / bases mapped (cigar)
average length: 6816
maximum length: 6816
average quality:    8.9
insert size average:    0.0
insert size standard deviation: 0.0
inward oriented pairs:  0
outward oriented pairs: 0
pairs with other orientation:   0
pairs on different chromosomes: 0

For a length filter, an awk length check can be slotted in instead of the grep:

$ samtools view -h mapped.bam |
    awk -F'\t' '{if((/^@/) || (length($10)>500)){print $0}}' |
    samtools stats | grep '^SN' | cut -f 2-

raw total sequences:    131
filtered sequences: 0
sequences:  131
is sorted:  1
1st fragments:  131
last fragments: 0
reads mapped:   131
reads mapped and paired:    0   # paired-end technology bit set + both mates mapped
reads unmapped: 0
reads properly paired:  0   # proper-pair bit set
reads paired:   0   # paired-end technology bit set
reads duplicated:   0   # PCR or optical duplicate bit set
reads MQ0:  0   # mapped and MQ=0
reads QC failed:    0
non-primary alignments: 0
total length:   569425  # ignores clipping
bases mapped:   569425  # ignores clipping
bases mapped (cigar):   453884  # more accurate
bases trimmed:  0
bases duplicated:   0
mismatches: 113755  # from NM fields
error rate: 2.506257e-01    # mismatches / bases mapped (cigar)
average length: 4346
maximum length: 16909
average quality:    9.1
insert size average:    0.0
insert size standard deviation: 0.0
inward oriented pairs:  0
outward oriented pairs: 0
pairs with other orientation:   0
pairs on different chromosomes: 0


Are there computational tools to extract features of DNA sequences? I am looking for tools to extract features from short DNA sequences. For example, entropy, complexity, GC-content, etc.
I have found the generateFeatures.py script from the PyFeat repo, but is there love a more widely used source code or a standard way to extract features from sequences biopython or similar?
Also, I think that since there are $2N$ encoded bits in a sequence of $N$ nucleotides, we have at most $2N$ independent features that could be extracted.

In addition, I am curios if there are any transformers models for DNA sequences.

If you want to do this with biopython, the SeqUtils package could be a solution.


Identify non-coding regions from a genome annotation I have this GTF file and I use the command below on a Linux machine to extract the coding regions of the genome:
awk '{if($3=="transcript" && $20=="\"protein_coding\";"){print $0}}' gencode.gtf

How I could do the inverse and keep only non coding regions?
I need non coding regions of human genome, but for asking my question here I referred to coding parts too.
After googling I see I actually need enhancers, promoters, or other regulatory elements of human genome to find cancer driver genes placed in these regions.
Likely non-coding regions of genome are here:
https://elifesciences.org/download/aHR0cHM6Ly9jZG4uZWxpZmVzY2llbmNlcy5vcmcvYXJ0aWNsZXMvMjE3NzgvZWxpZmUtMjE3Nzgtc3VwcDQtdjMueGxzeA==/elife-21778-supp4-v3.xlsx?_hash=KQi5jfO3kT2c4Qw44j4Rg6YAyCBQilYuWHVYXcRDuuo%3D

Getting the non coding regions of a protein coding transcript, sounds like you are looking for UTR.
UTR has its own feature in the gtf file. So you can do this:
$ awk -v FS="\t" '$3=="UTR"' gencode.gtf

If the gtf file is compressed use this instead:
$ zcat gencode.gtf.gz | awk -v FS="\t" '$3=="UTR"'

BTW: Why are you using such an old release of gencode? The current version is v29.
As @Wouter tells you, the non coding region of a genome is the complement of the coding regions. Coding regions have its own feature in the gtf file. You can get them with $ awk -v FS="\t" '$3=="CDS"' gencode.gtf. Reading the manual for bedtools complement is your task.


Frequency of specific viral sequence in .BAM or .fastq I was wondering if anyone had any experience 'counting' the frequency of a particular viral sequence in an individual's fastq sequence. So basically, counting the number of occurrences of a particular sequence given someone's fastq file.

The thing is, we must speak with the sequencing company regarding what format we'd like the files and must choose between .bam and .fastq - both would be extra $. Usually I'd choose .bam out of convenience but I was wondering, given that it has already been aligned the human genome I might not be able to proceed with this 'counting' process. So, my question is: if I were to use samtools - or any other program - to convert the .bam files to .fastq, is there any possibility that I could lose such information - e.g. viral seqs - that would have otherwise been available should I have just gotten the .fastq files directly?

Most read aligners will report unaligned reads as well, which presumably will include your viral sequences. I would ask them to formally confirm that the BAM files will contain unaligned reads before choosing that option.


How do I increase the minimum dot size in Seurat's DotPlot function? I want to use the DotPlot function from Seurat v3 to visualise the expression of some genes across clusters. However when the expression of a gene is zero or very low, the dot size is so small that it is not clearly visible when printed on paper. Please is there a possibility to increase the minimum dot size in the DotPlot function to make the dot sizes more visible when printed?

For example, I would like to have a minimum dot size set to be like

from this publication

Rather than:

from this other publication

Thank you in advance for your helpful hint.

Did you try to use DotPlot(..., scale.by = "size")?

In contrast to the default scale.by= "radius", this will link the area (==2*pi*r^2), not the radius, of the circles to the fraction of cells expressing the feature.
This corresponds much better to our perception of size and will make differences in low values easier to see.

Reading ?Seurat::DotPlot the scale.min parameter looked promising but looking at the code it seems to censor the data as well.

Since Seurat's plotting functionality is based on ggplot2 you can also adjust the color scale by simply adding scale_fill_viridis() etc. to the returned plot. This might also work for size. Try something like:

DotPlot(...) + scale_size(range = c(5, 10)) # will like warn about supplying the same scale twice

Apart from this, Seurat's plotting system is not very hackable and I find it much easier to extract the relevant data and plot them myself with ggplot2.

Yet another comment: Your plot with the strong differences looks much more convincing to me wrt. to the marker property of these genese than thee cited plot.


How to filter out cross alignments from a BED file? I have a BAM file:

@SQ SN:chr1 LN:248956422
@SQ SN:chrx LN:248956423
ST-E00110:348:HGVKKALXX:1:1201:5822:48670   323 chr1    9999    0   67H66M16H   chrx    1000    0   GATAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAAC  JJJJJJJJJJJJJJJJAJJJJJJJJJJJJFJJJJJJFJFJJJJJJFJJJJJJJJJJJA77FJFJJJ  NM:i:0  MD:Z:66 AS:i:66 XS:i:65 SA:Z:chr5,18606834,-,73S76M,34,0;   RG:Z:g1

There is a read aligned to chr1, and it's mate aligned to chrx.

I have a BED file:

chr1    0   100000  TestOnly

I would like to filter out everything that falls outside my BED region, that includes cross-alignments. In my example, although my read aligned to chr1 but it's mate is not. I don't want this read.

When I do:

  samtools view -L test.bed test.bam

the command gives me the read because it doesn't check cross-alignments.

My solution:

  samtools view -L test.bed test.bam | grep -v chrx

but this is very slow and clumsy. In my production pipeline I would have to do something like:

  samtools view -L test.bed test.bam | grep -v chrx | grep -v ... | grep -v ... | grep -v ... | grep -v ...

Q: Is there a better solution?

According to the SAM specification, the 3rd field of a SAM line (RNAME) is:

  RNAME: Reference sequence NAME of the alignment. If @SQ header lines
  are present, RNAME (if not ‘*’) must be present in one of the SQ-SN
  tag. An unmapped segment without coordinate has a ‘*’ at this field.
  However, an unmapped segment may also have an ordinary coordinate such
  that it can be placed at a desired position after sorting. If RNAME is
  ‘*’, no assumptions can be made about POS and CIGAR.

And the 7th field is (emphasis mine, missing "to" theirs):

  RNEXT: Reference sequence name of the primary alignment of the NEXT
  read in the template. For the last read, the next read is the first
  read in the template. If @SQ header lines are present, RNEXT (if not
  ‘*’ or ‘=’) must be present in one of the SQ-SN tag. This field is set
  as ‘*’ when the information is unavailable, and set as ‘=’ if RNEXT is
  identical RNAME. If not ‘=’ and the next read in the template has one
  primary mapping (see also bit 0x100 in FLAG), this field is identical
  to RNAME at the primary line of the next read. If RNEXT is ‘*’, no
  assumptions can be made on PNEXT and bit 0x20

So, you want to remove those lines whose 7th field isn't = and, just in case, those lines whose 7th field isn't = and isn't the same as the 3rd field. You can therefore use something like this:

samtools view -L test.bed test.bam | awk '$7=="=" || $3==$7

And, to save as a bam file again:

samtools view -L test.bed test.bam | awk '$7=="=" && $3==$7 |
    samtolls view -b > fixed.bam

On a separate note, there's very rarely a need to chain multiple grep commands like that. You can just use \| (or | with the -E or -P options) to separate them. Something like:

samtools view -L test.bed test.bam | grep -v 'chrx\|chr2\|chr10\|chrN'

Or

samtools view -L test.bed test.bam | grep -Ev 'chrx|chr2|chr10|chrN'


Passing custom distance matrix to T_Coffee I am trying to pass a custom distance matrix to T_Coffee, as per the docs, but can't seem to get it to work.

The documentation says to pass a matrix like this using BLAST format
https://tcoffee.readthedocs.io/en/latest/tcoffee_technical_documentation.html#blast-format-recommended

# BLAST_MATRIX FORMAT
# ALPHABET=AGCT
A G C T
A 0 1 2 3
G 0 2 3 4
C 1 1 2 3
...

I've tried with and without the leading # comments, and with and without the trailing ..., but no matter what I keep getting this error:

Parsing_error when reading blast_matrix doolittle:
A C
0

Here, my distance matrix file is titled doolittle and is the difference between the hydropathy indices of amino acids based on the doolittle scale * 10 to make them all integers. I've pasted the full matrix below. I am calling the command from linux 18.08 ubuntu terminal (although I've tried a combination of others.)

t_coffee -in=rhomboids.fasta -matrix doolittle

Any help would be appreciated. Thank you!

doolittle

# BLAST_MATRIX FORMAT
# ALPHABET=ACDEFGHIKLMNPQRSTVWY
A C D E F G H I K L M N P Q R S T V W Y
A 0 7 53 53 10 22 50 27 57 20 1 53 34 53 58 26 25 24 27 31
C 7 0 60 60 3 29 57 20 64 13 6 60 41 60 65 33 32 17 34 38
D 53 60 0 0 63 31 3 80 4 73 54 0 19 0 5 27 28 77 26 22
E 53 60 0 0 63 31 3 80 4 73 54 0 19 0 5 27 28 77 26 22
F 10 3 63 63 0 32 60 17 67 10 9 63 44 63 68 36 35 14 37 41
G 22 29 31 31 32 0 28 49 35 42 23 31 12 31 36 4 3 46 5 9
H 50 57 3 3 60 28 0 77 7 70 51 3 16 3 8 24 25 74 23 19
I 27 20 80 80 17 49 77 0 84 7 26 80 61 80 85 53 52 3 54 58
K 57 64 4 4 67 35 7 84 0 77 58 4 23 4 1 31 32 81 30 26
L 20 13 73 73 10 42 70 7 77 0 19 73 54 73 78 46 45 4 47 51
M 1 6 54 54 9 23 51 26 58 19 0 54 35 54 59 27 26 23 28 32
N 53 60 0 0 63 31 3 80 4 73 54 0 19 0 5 27 28 77 26 22
P 34 41 19 19 44 12 16 61 23 54 35 19 0 19 24 8 9 58 7 3
Q 53 60 0 0 63 31 3 80 4 73 54 0 19 0 5 27 28 77 26 22
R 58 65 5 5 68 36 8 85 1 78 59 5 24 5 0 32 33 82 31 27
S 26 33 27 27 36 4 24 53 31 46 27 27 8 27 32 0 1 50 1 5
T 25 32 28 28 35 3 25 52 32 45 26 28 9 28 33 1 0 49 2 6
V 24 17 77 77 14 46 74 3 81 4 23 77 58 77 82 50 49 0 51 55
W 27 34 26 26 37 5 23 54 30 47 28 26 7 26 31 1 2 51 0 4
Y 31 38 22 22 41 9 19 58 26 51 32 22 3 22 27 5 6 55 4 0


In silico PCR: any models that can quantitatively predict the initial content from the end amplified result? Let's say I've run the PCR and I now have a certain amount of the amplicon product. Is there any model that infers how much of it was there before the PCR started? I suspect it should be a probabilistic model (a distribution) parameters (shape) of which will strongly depend on the details of the PCR protocol and the characteristics of the amplicon...

If not, then may be there are forward models that predict the amount of the amplified result from the initial content?

The basic equation for PCR kinetics (Equation 1) states that the amount of amplicon after c cycles (Nc) is the starting concentration of the amplicon (N0) times the amplification efficiency (E) to the power c. The PCR efficiency in this equation is a number between 1 and 2 (2 indicates 100% efficiency).  NCBI

In this  article  on research gate they said  amplifx  was the software you were looking for to, calculate "(i.e. expected efficiency in PCR reaction)". The amplifx even mentions "Testing primers in in silico PCR".


DESeqDataSetFromTximport all(lengths > 0) is not TRUE I am suddenly running into an error when running the DESeqDataSetFromTximport

txi.rsem <- tximport(files, type = "rsem", txIn = FALSE, txOut = FALSE)
dds <- DESeqDataSetFromTximport(txi = txi.rsem,
                                colData = SampleFile,
                                design = ~ Compound)`

using counts and average transcript lengths from tximport
Error in DESeqDataSetFromTximport(txi = txi.rsem, colData = SampleFile,  :
  all(lengths > 0) is not TRUE

any suggestions what causes this error?

I found this same error as of today using this code:
pad.gene.data <-
  tximport(files = file.list.for.gene.analysis,
           type = "rsem",
           txIn = FALSE,
           txOut = FALSE)

dds.gene.level <-
  DESeqDataSetFromTximport(txi = pad.gene.data,
                           colData = pad.characteristic.data,
                           design = ~ condition)

I had done this same differential gene expression analysis using transcript counts aggregated to the gene level, on a mac about 5 weeks ago, but as of this past week switched to a windows and re-installed R 3.6.3 and whatever the current version of tximport and DESeq2 is.
I was able to replicate the results I had gotten on the mac 5 weeks ago removing all rows where all gene lengths are greater than 0:
pad.gene.data$abundance <-
  pad.gene.data$abundance[apply(pad.gene.data$length,
                                1,
                                function(row) all(row !=0 )),]

pad.gene.data$counts <-
  pad.gene.data$counts[apply(pad.gene.data$length,
                             1,
                             function(row) all(row !=0 )),]

pad.gene.data$length <-
  pad.gene.data$length[apply(pad.gene.data$length,
                             1,
                             function(row) all(row !=0 )),]

So first removing any gene where there is a length count of 0 from the abundance and counts matrix, then removing the genes from the length matrix so that they are all the same size. It was able to run that way.
I am definitely skeptical of why I would need to do that as it doesn't seem that it was advertised, if this is a new version, that we need to preprocess the data beforehand. I hope someone who might know more of the code can comment on this.


kraken2 to OTU table creation I did run kraken2 and I get two kinds of output
How do I generate OTU table from these outputs for phyloseq usage?

You need to run Bracken to the Kraken2 report output to estimate abundance.


Subset data frame based on ID I tried to subset the dataframe in R since there are empty space after each line in the sequence which shows 'NA' when i read the files into R ,as a result of which when i subset it just takes the sequence Protein_IDS and sequence line by line not as whole

Im looking for this Protein_ID in one column and respective sequence in other

Proteins_IDs                Sequence
  <chr>                       <chr>
1 fig|6666666.151743.peg.5284 MERQAVAVVFHDWRRAIQQRNKTAAADRGVDPVGVGHNSKIESLQPRSSVCRASLIQVPG
2 NA                          RLLNGEIFYSLAEAKVLIEAWRRHYNTVRPHSSLGYRPPAPETATPPYPASGSASLHLRP
3 NA                          DMAAMGLIH
4 NA                          NA
5 fig|6666666.151743.peg.5285 MIDRQTAQVRRWSFVLGSSIAAMAAASATPAWAQCSPDPTVADSVTNCTGTDSDGLRVTT
6 NA                          DGSMVEVAPSASVSGAGAAAITVDIAPVSGSYYARKATINVGGSVDGGTQSGISVLSGSV

So i have made another file with list of protein ID which i want to subset from this file..I tried removing empty space from the file but couldn;t as it the sequence should come in one line but instead there is line break

data files

I tried this SEQ <- inner_join(seq,protein_id) whose output is this its like each ID is associated with multiple sequence

fig|6666666.151743.peg.5044 MTKNTIHRTILVCASLIAMPAMAQGTTEAGQEGASLADIVVTAQRRSENLQNVPISVTAI
2 fig|6666666.151743.peg.5063 MTVEQALARMLEGTGLEAQRLSPTAYRLVQAAPPASRRDPVRPAPPLPPEGSVDIIVTGT
3 fig|6666666.151743.peg.5135 MLRVADGYVEEAADDAIVPSPAVPLSVPSAPVRFEAAPARYGARRGPNVPAILVIVLLHA

I want to advise you to use FASTA format, instead of obscure unknown formats. In R, there are many packages that can handle FASTA format and it is therefore more convenient. I will show you how to make a FASTA file from your data, and then I will show an example of how to import such a FASTA file and sub select your proteins of interest.

input <- scan("seq2.txt", what = character(), sep = "\t")
length(input)
new <- gsub("f", ">f", input[-1:-2])
write.table(new, "seq2.fasta", quote = F, row.names = F, col.names = F)

library(seqinr)
fasta.file <- read.fasta("seq2.fasta", seqtype = "AA", as.string = TRUE)
prot.id <- scan("protein_id.txt", what = character(), sep = "\t")
selection <- fasta.file[names(fasta.file) %in% prot.id]
length(selection)
[1] 87

head(selection, 2)
$`fig|6666666.151743.peg.5044`
[1] "MTKNTIHRTILVCASLIAMPAMAQGTTEAGQEGASLADIVVTAQRRSENLQNVPISVTAITSEQAIKLGITQPDVIAQMVPGLQLNHSAWGLTPFLRGVGTNSVAIGNEPAVALFVDDVYLPTGSAGVFQLNSISSIEVLKGPQGTLFGRNATGGVIHIRTRDPSFDPSASATASYGNYDTISGQLYGSTGLTDTVAINFALQGIDQRKGWGVNGTTGNEANSMSGWGGRAKLLWNMGEETSLLVAGSYDRRSSDQGQPGRSIAGTFTRDGYSADAAGLGFYDVMSNLDTVQKSTTRSLSGKLTHDFGSARLVLISAFSKTVVDPVVSDSDGSPSPFFQTREFAMGGRTFTQEVQLLSPEGSAVKWILGAFYLHDRSFMDATYEGQAFGTSKSHTDAVQRTNSIAGFGQVTTEILPRTNLTLGLRYTSDRRSLDAVARLLPATGAAVTFGPFEDAKTFNSVTGRASIDHHFSDDLMGYIAYNRGFKSGIYNFTGIVAGTTTLPPPVNPETLDAYTIGFKSEFLNRRLRVNAEAFLYEYSDLQVQIVLPTGTTLTNASKATIKGIDFEISAIPVRNLTVSASIGLLDGRYKDFRNAPSFFLLPPNAPTPYDPARCSGPNPYPAVQGSAPATQRSCDVSGNRLIQAPPFSSNLQLLYTLPTDLGDFDLSASWAHGGNYEYDVFSSPLTRQPRTDIVNGSVQWTAPGGTFNVRLWANNLTGEKYYTIATTGSTAGAKYAAAAPRTYGVTLGANF"
attr(,"name")
[1] "fig|6666666.151743.peg.5044"
attr(,"Annot")
[1] ">fig|6666666.151743.peg.5044"
attr(,"class")
[1] "SeqFastaAA"

$`fig|6666666.151743.peg.5063`
[1] "MTVEQALARMLEGTGLEAQRLSPTAYRLVQAAPPASRRDPVRPAPPLPPEGSVDIIVTGTKREEPLARLPISVAIVQAPGEQGGALLPGTEEVAQRLDGLTLTNLGPGRNRQFIRGVADSPFNGNSQSTVAILINDARATYNAPDPDLKLVDVERVELLKGPQGSLYGTGTLGGVYHIVTRPPDAGSFAAMAGMGASLATHGEMGGSAEAMVNLPVVADRLALRAVAYAVREPGWIDTAGRGDANRGELRGGRVALRALAGDWTVDLNGMAQMRDVADSQYVYARHALERPAQPAEPHDNDFALMALRVAGKLGGADLVYAGNWVDHDVSSTLDASDLSASAFGIAPPALYLEQRRNRLFDQEVRLSQTRDDGLSWLAGLSFVDVRNHLSGVVRGVRDETVGTGIQGVVEFSAFGDMRMPLTGALSLSVGGRLFHSRIEEEQGARHARIEASVARSGFTPSLALSWQAGAHDFLYARYAGALRPGGLSRTGAASAERFESDELQTVEMGWRHSASAWLTVNLDAFLTTWSHMQSDYLLADGLIGTRNVGDARIRGVEASVDWKPAAGWAISAGLVGQSALLTRTNEGLRLEDRHLPVVPDWKARFGVNRHFDLSGWSGSMGVNLAYLGPSRLSFDPGLDRRMGDYALANADLALRRKDWSIRLSGSNLLDCKADSFAFGNPFALARGTQYTPLRPRQWAIGFTKNW"
attr(,"name")
[1] "fig|6666666.151743.peg.5063"
attr(,"Annot")
[1] ">fig|6666666.151743.peg.5063"
attr(,"class")
[1] "SeqFastaAA"


Getting protein FASTA sequence based on keyword with python I would like to gather proteins FASTA sequence from Entrez with python 2.7. I am looking for any proteins that have the keywords: "terminase" and "large" in their name. So far I got this code:

from Bio import Entrez
Entrez.email = "example@example.org"

searchResultHandle = Entrez.esearch(db="protein", term="terminase large", retmax=1000)
searchResult = Entrez.read(searchResultHandle)
ids = searchResult["IdList"]

handle = Entrez.efetch(db="protein", id=ids, rettype="fasta", retmode="text")
record = handle.read()

out_handle = open('myfasta.fasta', 'w')
out_handle.write(record.rstrip('\n'))

However, it can get me several terminases from various organisms, while I need only terminase form bacteriophages (specificly Viruses [taxid 10239], host: bacteria). I've managed to get the nuccore accession ids from NCBI of the viruses I am intersted in, but I don't know how to combine those two informations. The id file looks like this:

NC_001341
NC_001447
NC_028834
NC_023556
...

Do I need to access every gb file of every ID and search for my desired protein in it?

Found what I was looking for. In:

searchResultHandle = Entrez.esearch(db="protein", term="terminase large", retmax=1000)

I've added:

searchterm = "(terminase large subunit AND viruses[Organism]) AND Caudovirales AND refseq[Filter]"
searchResultHandle = Entrez.esearch(db="protein", term=searchterm, retmax=6000)

which norrowed down my searches to the desired viruses. Granted it's not filtered by host, but by a taxonomy group, but it is enough for my work.

Thank you @Llopis for additional help


How to interpret the SMRT Link base modification algorithm output? I have got PacBio data from a Sequel machine at a coverage > 50X. I have run it in SMRT Link with their Base Modification Detection algorithm. I am now trying to figure out how to interpret the output and where to find the information about the modified bases. One of the outputs is a GFF3 file. The top of this file looks like this:

> ##gff-version 3

...

> I kinModCall  modified_base   1168    1168    31  -   .   coverage=48;context=TTAAACTAATTTCCAAAAATTTTGAACAGTTTTTAATACAA;IPDRatio=2.96
> I kinModCall  modified_base   1497    1497    32  -   .   coverage=50;context=CTTGAAAACACGAATATCTGTTTAGACAATAGTTTTGCCAA;IPDRatio=3.19
> I kinModCall  modified_base   1509    1509    33  -   .   coverage=48;context=CATGAAAAATTTCTTGAAAACACGAATATCTGTTTAGACAA;IPDRatio=2.84
> I kinModCall  modified_base   1529    1529    39  -   .   coverage=46;context=AAGAAAAGACCAAGAAAAACCATGAAAAATTTCTTGAAAAC;IPDRatio=3.84
> I kinModCall  modified_base   1904    1904    36  -   .   coverage=54;context=TAAAACCTGACAGAAAGAATTGGTCCCGTCATCTACTTTCT;IPDRatio=2.86
> I kinModCall  modified_base   2152    2152    31  -   .   coverage=59;context=TTAAAAATTAATATAAAATCAGGGAAATTTTTTTAAATTTT;IPDRatio=2.98
> I kinModCall  modified_base   2226    2226    36  -   .   coverage=63;context=TGGGCGGGGCGTCTCTGCCGGGAAAATTTTGCGTTTAAAAA;IPDRatio=2.59
> I kinModCall  modified_base   2228    2228    34  -   .   coverage=52;context=CGTGGGCGGGGCGTCTCTGCCGGGAAAATTTTGCGTTTAAA;IPDRatio=2.98
> I kinModCall  modified_base   2237    2237    33  -   .   coverage=60;context=GCACGGTTTCGTGGGCGGGGCGTCTCTGCCGGGAAAATTTT;IPDRatio=3.07
> I kinModCall  modified_base   2240    2240    44  -   .   coverage=63;context=GCGGCACGGTTTCGTGGGCGGGGCGTCTCTGCCGGGAAAAT;IPDRatio=2.89
> I kinModCall  modified_base   2311    2311    31  -   .   coverage=63;context=AAAATCGACGAAAATCGGTATAAAATCAAATAAAAATAGAA;IPDRatio=4.35
> I kinModCall  modified_base   2313    2313    39  -   .   coverage=62;context=GAAAAATCGACGAAAATCGGTATAAAATCAAATAAAAATAG;IPDRatio=4.00
> I kinModCall  modified_base   2314    2314    34  -   .   coverage=63;context=AGAAAAATCGACGAAAATCGGTATAAAATCAAATAAAAATA;IPDRatio=2.63
> I kinModCall  modified_base   2333    2333    35  -   .   coverage=62;context=AAAAAAAAGAGAAAAAATGAGAAAAATCGACGAAAATCGGT;IPDRatio=3.77
> I kinModCall  modified_base   2334    2334    35  -   .   coverage=63;context=CAAAAAAAAGAGAAAAAATGAGAAAAATCGACGAAAATCGG;IPDRatio=2.77

I expected something along the lines of 'chromosome I position 1497 methylcytosine', etc.

So it seems that the algorithm has detected modified bases at the following positions but it cannot call them, maybe because of a lack of coverage?

Update: the screenshot below shows there is no CSV file in the Data tab of the completed analysis.

It sounds like you're looking for the modifications.csv file, rather than the modifications.gff file. The .gff file is for viewing in SMRT View (or other viewers). You can see schema for both files at PacBio's Methylome Analysis Technical Note. Specifically, modified_base means the algorithm could not determine the exact base and modification.

I also found a video tutorial for performing the analysis in SMRTLink on YouTube

EDIT: Per PacBio's Methylome Analysis Technical Note documentation, the .csv "...is also produced when motif analysis is not active, such as with the RS_Modification_Detection analysis protocol." Not sure why it's missing, perhaps because the algorithm couldn't clearly identify any of the modifications. You'll probably need to contact tech support on this issue.

In either case, the .gff you posted includes the significant hits (p < 0.01), but the algorithm was unable to classify them, likely because of insufficient coverage. PacBio's Methylome Analysis Technical Note documentation gives specific coverage recommendations for each event type, but 100x is the lower bound. Lots of good information there.

Probably not what you wanted to hear, but I hope that gets you going.


Python module for fetching NCBI id for a list of species I have a list of scientific names of species. Is there a python module that can fetch NCBI taxonomy IDs?

Biopython has an efetch frontend.

from Bio import Entrez
Entrez.email = "michaelg@example.ac.uk"
handle = Entrez.efetch(db="nucleotide", id="KR052012.1", rettype="gb", retmode="text")
print(handle.read())

OUT

LOCUS       KR052012               10714 bp    RNA     linear   VRL 09-MAY-2015
DEFINITION  Dengue virus strain Hb33/CHN/2014, complete genome.
ACCESSION   KR052012
VERSION     KR052012.1
KEYWORDS
.... etc

You will likely need to perform an EFetch search and then extract out the accession numbers for import (as done above)


filter secondary alignments using pysam I am new to python and trying to learn. The below is an attempt to filter out secondary reads in a bam using pysam. I ran mark duplicates and got an error on several bam files due to secondary alignments. I added comments to each line to illustrate my thought process.  Thank you :).

The specific mark duplicates error:

Value was put into PairInfoMap more than once

Here is the script:

#! /usr/bin/python  ## call python script
import sys ## import python system functions
import pysam  ## import module

bam = pysam.AlignmentFile(*.bam, "rb")  ## open bam and read
click.echo("Reading BAM file")  ## output message
for read in bam:  ## start loop and iterate over each bam
    if read.is_secondary=true  ## not the primary alignment
    if read.has_tag('XS')  ## look for Secondary Alignment Score tag
        XS = read.get_tag('XS')  ## read and store AS value
    if read.has_tag('AS')  ## look for Alignment Score tag
        AS = read.get_tag('AS')  ## read and store value os XS tag
    if AS < XS  ## Alignment score less then
        read.write(read)  ## only print primary alignments
bam.close()  ## end processing

Quick look at your script suggests that you are comparing the XS and AS from different reads. I would recommend removing the continue from the read.has_tag('XS') continue block.

I think that if you filter out is_secondary and is_supplementary you don't need to do the XS to AS comparison assuming whatever generated for BAM did so correctly.

for read in bam.fetch()
    if read.is_secondary=true
        continue
    if read.has_tag('XS') <- This block stores XS and then moves on to the next read
        XS = read.get_tag('XS')
        continue
    if read.has_tag('AS')
        AS = read.get_tag('AS')
    if AS < XS  <- This block compares the AS of the current read to the XS of the previous read that had one.
        read.write(read)
bam.close()


How to find all WGS assemblies accessions of a species Some background
Similar to the OP of https://www.biostars.org/p/377840/, I would like to programmatically BLAST a sequence to a local database of all WGS assemblies.
Since this isn't feasible for the average biology lab server (correct me if I am wrong), I plan to use ncbi-acc-download to download all WGS assemblies of the species of interest (not a popular species like E. coli, so it should be feasible). Then, I will create a BLAST database for the downloaded assemblies and BLAST the sequence to it.
My question
How can I find all WGS assemblies accessions of a species?

My current plan is to search the NCBI Assembly database using Entrez and a search term such as "wgs"[Properties] AND txid1337[orgn:exp].
EDIT: IIUC, this approach might miss some WGS assemblies. See my answer.
I am worried (and thus ask for your help) this isn't the right approach because there seem to be at least 3 other places in which assemblies can be found:

The assembly_summary.txt files that are described in https://ftp.ncbi.nlm.nih.gov/genomes/README_assembly_summary.txt.
https://www.ncbi.nlm.nih.gov/Traces/wgs/. This is the one I am most worried about. Both because it seems to have much more entries, and also because https://www.ncbi.nlm.nih.gov/genbank/wgs/ says: "The public WGS projects are at the list of WGS projects ."
taxid2wgs.pl. (As mentioned in https://www.biostars.org/p/377840/#471535, this one seems to miss some assemblies.)

There seem to be WGS assemblies that can't be found in NCBI Assembly database, e.g.: https://www.ncbi.nlm.nih.gov/nuccore/1779902990. I guess that such assemblies also cannot be found in The assembly_summary.txt files that are described in https://ftp.ncbi.nlm.nih.gov/genomes/README_assembly_summary.txt.
My current best guess is that each WGS assembly has a "WGS master record" in NCBI Nuccore database. To find all WGS assemblies of a taxon whose uid in NCBI Taxonomy database is 1337, search NCBI Nuccore database using Entrez and the search term "wgs master"[Properties] AND txid1337[orgn:exp].


Deep learning RNA sequences Currently I'm working on a project, which combines deep learning with RNA sequences. I'll try to predict pseudotorsion angles [1] from raw rna sequence. The ideas is to train a neural network with raw rna sequences and for each nucleotide their corresponding pseudotorsion angles, and than to predict the angles of the remaining sequences in the test set.

How my data is structured:
example:

seq1:   A  C  G    G   U   A   C
Eta:   169 87 110 87  45  187 78
Theta: 123 10 45  168 132 34  100

[1] These are angles describing the backbone conformation of a rna molecule.

I'm pretty new to the field of deep learning, and so far I build a simple feedforward neural network, but it's prediction accuracy is pretty low with only one percent.

Has anyone some tips for me how to improve this?
How do I preprocess this kind of data correctly for deep learning?

I appreciate any help.

The question of folding of RNA sequences looks slightly similar to protein folding - perhaps searching in this domain might bring more suggestions.

An example of (current) state of the art of deep learning approach to protein folding is AlphaFold developed by DeepMind for CASP competition . They basically decompose the angles to an image representation and look for structures in the proximity of each atom. While it requires a significant processing power, RNA structures might be easier due to the size of the sequence.


Error: FASTA file not found using building kallisto index I am trying to build an index on kallisto on mac using R studio and have saved my FASTA file to my desktop, however, when I run kallisto index -i DI_mRNA.index DI_mRNA.fa.gz, there is an error saying Error: FASTA file not found DI_mRNA.fa.gz

As mentioned in a comment, your working directory lacks the fasta file.


Database for germline copy number variations in cancer I am interested in looking at germline copy number variations in individuals that are at high risk of developing cancer. Are there any databases where I can look, if a CNV seen in our test case has been seen/reported earlier in a patient with cancer?
Thanks.

You can find CNV (and other structural variant) information in:

DECIPHER

  DECIPHER is used by the clinical community to share and compare phenotypic and genotypic data. The DECIPHER database contains data from 29603 patients who have given consent for broad data-sharing; DECIPHER also supports more limited sharing via consortia. Have a look at the numbers.

  Anyone can browse publicly-available patient data on DECIPHER and request to be put in contact with the responsible clinician. Why? Because sharing benefits everyone.

ExAC

  The Exome Aggregation Consortium (ExAC) is a coalition of investigators seeking to aggregate and harmonize exome sequencing data from a wide variety of large-scale sequencing projects, and to make summary data available for the wider scientific community.

  The data set provided on this website spans 60,706 unrelated individuals sequenced as part of various disease-specific and population genetic studies. The ExAC Principal Investigators and groups that have contributed data to the current release are listed here.

ExAC has included CNV calls since August 2016.
ClinVar

  ClinVar is a freely accessible, public archive of reports of the relationships among human variations and phenotypes, with supporting evidence. ClinVar thus facilitates access to and communication about the relationships asserted between human variation and observed health status, and the history of that interpretation. ClinVar processes submissions reporting variants found in patient samples, assertions made regarding their clinical significance, information about the submitter, and other supporting data.

Database of Genomic Variants

  The objective of the Database of Genomic Variants is to provide a comprehensive summary of structural variation in the human genome. We define structural variation as genomic alterations that involve segments of DNA that are larger than 50bp. The content of the database is only representing structural variation identified in healthy control samples.

  The Database of Genomic Variants provides a useful catalog of control data for studies aiming to correlate genomic variation with phenotypic data. The database is continuously updated with new data from peer reviewed research studies. We always welcome suggestions and comments regarding the database from the research community.

Information from all of these is also displayed in the Structural Variant browser of VarSome*. While you can't search by patient status, so you can't use VarSome to know the clinical history of any sample with the CNV, you can at least look up a specific gene and see the structural variants listed for it. For example, for TP53, you would see:

* Disclaimer: I work for the company behind VarSome.


How to input data for DESeq2 from individual HTSeq count? I am comparing the gene expression of 2 bacteria under 1 condition. I have now the count tables for 3 tech. replicates for each bacteria.

Bacteria1_1.count
Bacteria1_2.count
Bacteria1_3.count

...same for the other bacteria.

These files look like this:

gene1 10000
gene2 500
gene3 0
gene4 5000

I want to use DESeq2 for differential gene expression analysis. But I cannot figure out how to properly execute the DESeqDataSetFromHTSeqCount() command with this type of data.

Is there another intermediate step to add?

In case still needed, here is a piece of code that I have used for the same purpose (includes a lot of sanity checks). This is highly likely an overkill as each count file would probably have the same rows (as a result of being generated from the same workflow) and thus can easily be merged column-wise.

count_file_dir <- "../input/count_files/"
count_file_names <- list.files(count_file_dir)
count_files <- paste0(count_file_dir, count_file_names)

# read and merge count files, while removing last 5 lines that correspond to summary stats
count_table_list <- mclapply(count_files, function(x) {fread(x)[1:(.N-5)]})
names(count_table_list) <- sub(".count", "", count_file_names)

# check the dimensions of count files
t(sapply(count_table_list, dim))

# check if the count files have the same gene ordering
# compare gene names of the first file with all
gene_names <- sapply(count_table_list, function(x) x$V1)
as.data.frame(apply(gene_names, 2, function(x) all(x == gene_names[,1])))

# add file name as a unique column name for gene counts
for(i in names(count_table_list)) {
  names(count_table_list[[i]])[2] <- i
}

# merge read counts and convert to a matrix
count_matrix <- Reduce(merge, count_table_list)
count_matrix_rownames <- count_matrix$V1
count_matrix$V1 <- NULL
count_matrix <- as.matrix(count_matrix)
rownames(count_matrix) <- count_matrix_rownames
colnames(count_matrix) <- sub(".count", "", colnames(count_matrix))

dim(count_matrix)

dds <- DESeqDataSetFromMatrix(countData = count_matrix,
                              colData = ...,
                              design = ...)


Input from multiple channel, merge and use in nextflow I need to merge multiple channel to reuse them in two processes. So far I tried as follows:
channel.from([15,20,21]).into {k1;k2}
channel.from(['a','b','c']).into {m1;m2}
channel.from([1,2,3]).into {l1;l2}
para1=m1.merge(k1).merge(l1)
para2=m2.merge(k2).merge(l2)
process p1{
input:
each m, k, l from para1
"""
echo $m, $k, $l
"""
}
process p2{
input:
each m, k, l from para2
"""
echo $m, $k, $l
"""
}

The error I am getting:
Unknown process directive: '_in_each'

Did you mean of these?
    each

If I use para1.view() I can see the merged elements as expected. Any help

I don't think you can unpack tuples using the each qualifier like that. You can, however, unpack them in the script block:
process test {

    input:
    each triple from para1

    script:
    def (m, k, l) = triple

    """
    echo "${m} ${k} ${l}"
    """
}

Consider also that the merge operator is deprecated and will be removed from future Nextflow releases. Instead, you may need to join by index, for example:
def indexedChannel( items ) {
    return Channel.from( items.withIndex() ).map { item, idx -> tuple( idx, item ) }
}

k_ch = indexedChannel( [ 15, 20, 21 ] )
m_ch = indexedChannel( [ 'a', 'b', 'c' ] )
l_ch = indexedChannel( [ 1, 2, 3 ] )

m_ch
    .join( k_ch )
    .join( l_ch )
    .map { idx, m, k, l -> tuple( m, k, l ) }
    .into { para1; para2 }


How to fit all genes (labels) in chromosome ideogram plot made by RCircos package? I am using Rcircos to make a chromosome ideogram plot for my gene list (n = 45). However, I am getting this error:
 Not all labels will be plotted.
 Type RCircos.Get.Gene.Name.Plot.Parameters()
 to see the number of labels for each chromosome.

I am using the below-given script and what to know is there any way to fit all my genes in the chromosome ideogram plot
 library(RCircos)
 data(UCSC.HG19.Human.CytoBandIdeogram)
 RCircos.Set.Core.Components(
 cyto.info=UCSC.HG19.Human.CytoBandIdeogram,
 chr.exclude = NULL,
 tracks.inside = 10, tracks.outside = 0)

 #Making a plot with RCircos
 out.file <- "RCircos_IF.pdf";
 pdf(file=out.file, height=8, width=8, compress=TRUE);
 RCircos.Set.Plot.Area();
 par(mai=c(0.25, 0.25, 0.25, 0.25));
 plot.new();
 plot.window(c(-2.5,2.5), c(-2.5, 2.5));

 # plot chromosome ideogram
 RCircos.Set.Plot.Area();
 RCircos.Chromosome.Ideogram.Plot()
 gene.lable.data = read.table ("IF.txt", sep="\t", header=T, row.names = 1) # IF.txt is my input file
 RCircos.Gene.Connector.Plot(gene.lable.data, track.num = 1, side = "in");
 track.num <- 2;
 RCircos.Gene.Name.Plot(gene.lable.data,name.col = 4, track.num)

 dev.off()

My input file looks like
    Chromosome  chromStart  chromEnd    Gene.name
    1   chr1    150363091   150476566   RPRD2
    2   chr1    150549369   150560937   ADAMTSL4
    3   chr1    91949371    92014426    BRDT
    4   chr1    31365625    31376850    FABP3
    5   chr1    150960583   150975004   CERS2

To change the maximum number of genes in RCircos, we could modified the char.width according to this link.
The actual value in a RCircos session could be modified with get and reset methods for plot parameters.
params <- RCircos.Get.Plot.Parameters() #$char.width
params$char.width <- 100 #default 500
RCircos.Reset.Plot.Parameters(params)

#the maxLabels are updated accordingly
RCircos.Get.Gene.Name.Plot.Parameters()


How do I pull singe cell RNA sequencing data from GEO database? I am new to R and computational biology. I am trying to look through a published data set to check for gene expression for my own project. I am having trouble finding materials to teach me how to analyze single cell RNA sequencing data on R from the GEO NCBI database.

This is the data that I am trying to access https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE138651

I would appreciate any directions or suggestions for tutorials.

Thanks!

The simplest would be using a count matrix (at the end of the link you have shared, section "Supplementary file"). For example GSE138651_Vagal_WholeNodoseSeq_raw.mtx.gz corresponds to "raw" count matrices obtained after mapping and counting the reads. This file, along with the GSE138651_barcodes.tsv.gz and GSE138651_genes.tsv.gz can be used to "read" the data into R; to that end there are many options in the form of packages, Seurat is one of the best documented and hence is beginner-friendly.


Pseudo-temporal ordering in heterogeneous populations What's the most widely accepted tool for doing pseudo-temporal ordering from scRNAseq data? Also is there away to separate differential expression that occurs based on "cell identity" or maybe more accurately cell type fate from that which arises from cells being in different stages differentiation.

To be more concrete, lets say there's a population of cells, some of which were born at time 1, time 2, and time 3. The progression along the temporal trajectory can be described via a set of genes that are fluctuating as the cell matures. So you might have the same cell type which was born at time 1 be transcriptionally distinct from a younger one born at time 3.  On the other hand, within this population there are subpopulations which will have different cell fates and are transcriptionally distinct. Is there away to reliably separate the temporal axis from the "cell fate axes". If not is this something people are working on or is flawed logic to think this kind of thing is possible?

Regarding the first part of the question: scRNA-seq is a rapidly developing field so it may be hard to talk about "widely accepted tool for doing pseudo-temporal ordering from scRNAseq data". Few of the tools aiming to do just that include Monocole, Waterfall or Sincell (see this paper for references https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4122333/)

The second part is a bit more complex. Many aspects have come together here, e.g. cell-cycles, sub-populaitons identification and pseudo-temporal ordering, to get a true reflection of the biological processes. There are efforts on all (and more) these fronts (see again e.g. the above paper, though not the latest) and there are most likely people working on their integration. I'm not aware of any studies published yet at this depth though


alpha diversity wilcox.test I am trying to do Wilcox test to detect the significant difference in alpha diversity but it is showing an error? physeqN2 is a phyloseq object and Season is a metadata column.
results = estimate_richness(physeqN2, measures = 'Chao1')

d = sample_data(physeqN2)

# calculate Wilcox Test
Wet.root = results[d[,'Season'] == 'Wet',]
Dry.root = results[d[,'Season'] == 'Dry',]
pmv <-pairwise.wilcox.test(Wet.root, Dry.root)

capture.output(pmv, file = "chao1.16s.wet_vs_Dry.root.txt")

Error;
Error in wilcox.test.default(xi, xj, paired = paired, ...) :
  'x' must be numeric
Calls: pairwise.wilcox.test ... compare.levels -> wilcox.test -> wilcox.test.def                                                                                ault
Execution halted

I hvae tried to do it in another way but here also getting error?
alpha.diversity <- estimate_richness(physeqN2, measures=c("Chao1"))

## Wilcox on observed richness
data <- cbind(sample_data(physeqN2), alpha.diversity)
pmv <-pairwise.wilcox.test(Chao1 ~ Season, data)

capture.output(pmv, file = "chao1.16s.wet_vs_Dry.root.txt")

Error:

Error in if (length(ans) == 0L || as.character(ans[[1L]])[1L] == "~") { :
  missing value where TRUE/FALSE needed
Calls: pairwise.wilcox.test ... sapply -> lapply -> FUN -> compare.levels -> [ -> [.formula
Execution halted

Many thanks

I think this is a bug in the input data.
At a wild guess what I think is happening is the input data is not regarded as paired, hence the error paired = paired. The other weird error is x must be numeric. Firstly,

how are you scoring missing data, its not via - by any chance? Missing data must be via nan type syntax.
are you sure there are no heading bugs? Like a header is getting dragged into the numeric data?

Wilcoxan is a non-parametric T test, so the wet root and dry root must match exactly the same sample. Thus it is not a general wet root versus generic dry root, but wet and dry from the same specimen. If this is the data you have the program doesn't think it is matching.
If it is truely paired, again I think the source of error is you have some missing value which is leaving the pairings unequal lengths or else a header bug.
If the data genuinely isn't paired, then I would look at the Mann-Whitney U test, it would be easier in this scenario.
I don't use R syntax (just pandas dataframes) so can't say whether the syntax in the first example is valid, but it is pretty clear regarding your input. In Python you'd generally just give it the entire dataframe and tell it to go figure.
The reality is that loads of programs do Wilcoxan, you don't need to use R if its a one-off test. Obviously if this a part of a pipeline (automation) thats a bit different.


Using cellranger mkgtf w/ NCBI GFF3 files I am using several cellranger tools/pipelines to analyze some data, specifically, mkgtf, mkref, count, and aggr. I am using mkgtf and mkref to include a PI-synthesized protein to my reference. I have finished the analysis but am not seeing the protein in my output, so I am going back to the start to verify the process.

Unlike the documentation (https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/advanced/references), my GTF file was not downloaded from ENSEMBL or UCSC, but NCBI (https://www.ncbi.nlm.nih.gov/nuccore/AY678269.1?report=genbank). This GFF file does not have the same key nomenclature as seen in ENSEMBL or USCS.

I did use the mkgtf tool, but the created file is identical to the input.

cellranger mkgtf custom_ref/AY678269_genes.gtf reference_data/AY678269_mkgtf_genes.gtf

I went ahead and used cat to add the fasta and gtf file from the synthesized protein to the respective reference fasta and gene files, then ran mkref.

The s.protein is in the fasta file:

>lcl|AY678269.1_cds_AAV52169.1_1 [protein=tandem-dimer red fluorescent protein] [protein_id=AAV52169.1] [location=1..1431] [gbkey=CDS]

as well as the gene list

##sequence-region AY678269.1 1 1431
##species https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=32630
AY678269.1      Genbank region  1       1431    .       +       .
AY678269.1      Genbank CDS     1       1431    .       +       0

However this protein is not found in the output...surely I must have made an error. Does anyone have any ideas about what I could be doing incorrectly?

Do I need to change the entire nomenclature of the GFF keys to match what I see from ENSEMBL (e.g., gene_id ""; gene_version ""; transcript_id, etc...)?

Any assistance is greatly appreciated.

Regarding your comment:

  This GFF file does not have the same key nomenclature as seen in
  ENSEMBL or USCS

The link you have provided corresponds to an entry in Genbank format and not to GFF format. I believe you need a file in GTF format in particular.

Nevertheless, given the fact that you have manually added the corresponding entries to the reference fasta and gtf files prior to mkref, you should have seen your gene of interest in the output. I would double check one of the original fasta headers and original gtf entries in the reference file and match these for the new entry.


Can't align 2 sequences, (MAFFT killed) I have 2 .fa sequenes like:
>some description of first...
ACTGACTACGTACGACTACGATCTGACTACACGTAGATAGACTAGTCACTACG
ACTGACTACGTACGACTGACTACATGACTGAC...
> some description of second...
ATCGATAGCTACGATCGATAGTAGATAG....

This file is 461mb (1 chromosome of homo sapiens & 1 chromosome of pan troglodytes) and I want to align it. I am using mafft.
mafft --globalpair --maxiterate 1000 Homo_sapiens_1_Pan_troglodytes_1.fa > ./mafft_results/Homo_sapiens_1_Pan_troglodytes_1_output.fa

And it returns me
generating a scoring matrix for nucleotide (dist=200) ... done
All-to-all alignment.
/usr/bin/mafft: line 2719:  6030 Killed                  "$prefix/tbfast" _ -u $unalignlevel $localparam -C $numthreads $seqtype $model -g $pgexp -f $pggop -Q $spfactor -h $pgaof -A $usenaivepairscore $focusarg _ -+ $iterate -W $minimumweight -V "-"$gopdist -s $unalignlevel $legacygapopt $mergearg $outnum $addarg $add2ndhalfarg -C $numthreadstb $rnaopt $weightopt $treeinopt $treeoutopt $distoutopt $seqtype $model -f "-"$gop -Q $spfactor -h $aof $param_fft $localparam $algopt $treealg $scoreoutarg $focusarg < infile > /dev/null 2>> "$progressfile"

Why? How could I align it for further comparing?
EDIT 1
MAFFT working cool, when I tried to use it on test sequence file
>test1
ATCGATCGATGACTACGATCGACTACGTACG
ACTAGCTAGCTACGCATAGTAGATCGA
ATCGACTAGCTCAGATCGACTAGATCGATCGT
>test2
TTATCCACTGATCGACTCAGATCGACTCAGAC
ATCGATCAGCTACGATGACTACGTAC
ACTGACTAGCTAGACTGACTAGCTACGTCA

and after MAFFT mafft test.fa > output.fa it returned

test1
------atcgatcgat-------------gactacgatcgactacgtacgactagctagc
tacgcatagtagatcgaatcgactagctcagatcgactagatcgatcgt
test2
ttatccactgatcgactcagatcgactcagacatcgatcagctac---------gatgac
tacgtac----------actgactagct-agactgactagctacgtca-

So it works.
I can use all 4 at the example. Now about my large file:
alex@alex:~/Desktop/test/BOTH/homo_sapiens_pan_troglodytes$ mafft --thread 4 --maxiterate 2 Homo_sapiens_1_Pan_troglodytes_1.fa > ./mafft_results/Homo_sapiens_1_Pan_troglodytes_1_output.fa
nthread = 4
nthreadpair = 4
nthreadtb = 4
ppenalty_ex = 0
reallocating...
done.
reallocating...
done.
reallocating...
    done.
reallocating...
done.
reallocating...
...
stacksize: 8192 kb
reallocating...
done.
reallocating...
done.
reallocating...
done.
...
generating a scoring matrix for nucleotide (dist=200) ... done
Gap Penalty = -1.53, +0.00, +0.00

Making a distance matrix ..
    1 / 2 (thread    1)
done.

Constructing a UPGMA tree (efffree=1) ...
    0 / 2
done.

Progressive alignment 1/1...

len1=230480994, len2=223796828, Switching to the memsave mode
STEP     1 / 1 (thread    0)m/usr/bin/mafft: line 2719:  3192 Killed                  "$prefix/disttbfast" -q $npickup -E $cycledisttbfast -V "-"$gopdist -s $unalignlevel $legacygapopt $mergearg -W $tuplesize $termgapopt $outnum $addarg $add2ndhalfarg -C $numthreads-$numthreadstb $memopt $weightopt $treeinopt $treeoutopt $distoutopt $seqtype $model -g $gexp -f "-"$gop -Q $spfactor -h $aof $param_fft $algopt $treealg $scoreoutarg $anchoropt < infile > pre 2>> "$progressfile"

Doesnt matter I do not enter --thread 4 or not, --maxiterate 2 or --maxiterate 1000, it still shows me same error. It takes around 2 minute. Ram usage when running at the first minute is 6656mb, after 1 min 11028mb and then it out of ram and error..
PS
discussion moved here

--thread 4 is essential, --maxiterate just try and get it working, not getting it to work optimally. Just make sure there are no other processes running on the machine when --thread 4 is operating, e.g. you are using another machine as your desktop
The gist is as you spread the calculation across your processors (4) and simplify the calculation (lowering the iterations), you are getting further into the calculation, for example now you are getting into gap penalties.
If you are agreed then I would try performing the alignment using just one chromosome of the human genome against the chimp genome and I suspect that would work. Then take it from there.
The objective is to get the calculation working, for example I wouldn't trust a result with max. iterations of 2, but once you have it working you can then start refining the parameters towards optimal results. The i3 processor isn't perfect in context.
Using --thread 4 is super important and you want to increase maxiterate towards the default if the machine supports it.
I suspect the solution is to perform the calculation on an 8+-core machine and set --thread 8 (or higher depending on the number of cores on the machine), at least it will get further into the calculation. With 4-cores and an i3 (diagnostics of RAM usage would be nice) in any case it just doesn't seem to like it. The alternative is to perform this across the genome via sampling the chromosomes on your current i3.

Just to explain why this works. When you have a RAM intensive job and you spread it across cores the RAM overhead massively reduces for each core that is added. Don't ask me how this works in the hardware mecahanics of RAM distribution per core. The reason I know this is because of emperically parallelising Bayesian calculation using each chain on a different core.


GDC data portal cancer subtype I need to use Gene Expression quantification data (using HTSeq-FPKM-UQ) for breast cancer data on the gdc portal (https://portal.gdc.cancer.gov/repository) to classify breast cancer subtypes (luminal A, luminal B, HER2-like, basal-like).

However, since I don't know the portal, I'm having problems in associate the EXP data and the label.

I could get the label either by using the old TCGA-barcode and use a dictionary I already have from an article (but I don't know how to get the barcode from the portal / API) or by using clinical data (but I don't know how to find metadata that link the EXP and the clinical data).

Is there anyone more used to the portal who could help me?

I personally find the GDC portal not very user-friendly. However, some of the data is available from other sources. For example, you can download TCGA-BRCA FPKMs and phenotypes from UCSC Xena here: https://xenabrowser.net/datapages/?cohort=GDC%20TCGA%20Breast%20Cancer%20(BRCA)

All the sample info is already combined and the labels are the classic TCGA barcodes.


design formula question I am not sure I am building the proper design formula for the question I want to answer

I have the following samples with three factors; clone, the structure and the condition.

clone structure diabetic
1 07        2D      Dia
2 21        2D     Ctrl
3 23        2D      Dia
4 32        2D     Ctrl
5 34        2D      Dia
6 43        2D     Ctrl
7 07        3D      Dia
8 21        3D     Ctrl
9 23        3D      Dia
10 32        3D     Ctrl
11 34        3D      Dia
12 43        3D     Ctrl

I want to pull the differentially expressed genes for structure (2D vs 3D) blocking by clone, and the one differentially expressed by condition
 (regardless of the structure variable).

I am not sure how to build the design in one call, therefore I have produced the first object to pull 2D vs 3D blocking by clone using the following.

des_structure <- DESeqDataSetFromHTSeqCount(sampleTable = samples,
                                       directory = "../data/htseq_geneCounts/",
                                       design= ~ clone+structure)

for the condition controlled by structure I have used

des_condidtion <- ddsHTSeq_diabetic <- DESeqDataSetFromHTSeqCount(sampleTable = samples,
                                                 directory = "../data/htseq_geneCounts/",
                                                 design= ~ structure+diabetic)

I am not sure whether this is the best way to build the designs.

You can't (easily) use a single design, since clone is nested within diabetic status, so you were correct using two separate designs. Your designs are correct.

If you really wanted to use a single design for some reason, it'd be ~clone+structure and for the latter comparison you'd end up using a rather complicated contrast (something like c(rep(c(1,-1), 3), 0)). I really can't recommend that unless you have a very compelling reason, since it's vastly easier to think in terms of factors for your experiment and you're not hurting anything by fitting twice.


Single Cell RNA seq Analysis for low input Cells? I am having a single cell RNA seq data from ~500 Cells. I do understand this number is very low for present day analysis tools like Seurat. So, I want to know what if one have a such a small number of cells which tool works best?

Thanks

You can perfectly use any of the single cell processing packages out there. Seurat, Monocle3, Scanpy, Kallisto, etc. Following some of the tutorials and vignettes will most likely help you figuring out how to proceed.

These packages provide useful functions to explore single cell data, process the data and each of them has a series of advanced features to address specific questions. You don't need 10k cells to look at number of gene/cell, transcript/cells, % mitochondrial genes, etc. You can perfectly normalize and standardize 500 cells and run PCA, and analyze the complexity that describes your data.

To address if 500 cells is enough for you pilot we need to know what question you want to address.
The following tool for example gives some intuition on how many cells to sample if you wish to detect n-cells per cluster.

https://satijalab.org/howmanycells

Hopes this helps


Using Entrez.efetch() to retrive .fasta file from any NCBI database? I tried to use these python codes to retrieve a .fasta file from gene database, but the output were the parameters of the report
from Bio import Entrez
Entrez.email = "send_helps@earth.com"
handle = Entrez.efetch(db = "gene",id = "682",rettype="fasta", retmode="xml")
print(Entrez.read(handle))

i also read the documentation on this site about efetch() and seems like this function can only retrieve certain kinds of data for each database. So i can't take fasta file from gene database !?
I don't want to use the NCBI web UI. What's wrong with my code? Is there a way i can guarantee to retreive a .fasta file through python from NCBI ? Am i misundertanding something severly ?
Sorry for bad English. thank you.

In order to download the FASTA sequence of transcript variants using Entrez eutilities, you will have to use Entrez.elink to get a list of transcript accessions for the gene of interest and then use Entrez.efetch to fetch the sequence(s). Just to give you an idea, you can use Entrez Direct for this as follows:
elink -db gene -id 682 -target nuccore -name gene_nuccore_refseqrna \
  | efetch -format fasta

You can find more Entrez Link descriptions here. Note, there is a gene_fasta option for efetch but it still returns only the transcript sequence (i.e., no intron sequence).
If you need the full gene sequence (including the intronic regions), you can use the commands mentioned here.
Alternatively, you can use NCBI Datasets for this. You can use the web interface, command line tool or the python API. This page with the how-to articles has the information you seek.


Run Rscript for multiple samples (files) I have used the following script for calculating the median for one file (gene). However, I have 20,000 files for which I want to run this R script and get 200 separate output files

 #!/usr/bin/R
 setwd("/location folder")

 d <- read.table("tgf.txt", header=FALSE)

 x <- aggregate(d[, 2], list(d$V3), median)

 write.table(x, "tgf.output.txt", quote=FALSE, row.names=TRUE)

 q()

Please do not use setwd in a script, it changes working folder without explicit specifications from the user, so the results could be unpredictable for people that are not familiar with the content of the script.

Create a tab-separated file with input filename and output filename, or at least a file with all input files, unless you can pick them using a list.files(pattern = "<pattern_here>"). Read this file as a data.frame and use it with an apply, so a function processes the df row-by-row. Let's say the file (named metadata.txt) is:

in_file out_file
tgf.txt tgf.output.txt
abc.txt abc.filtered.txt
xyz.txt xyz.randomname.txt

Define a function called process_row() in the Rscript:

process_row <- function(r) {
    d <- read.table(r[[1]], header = FALSE);
    x <- aggregate(d[, 2], list(d$V3), median);
    write.table(x, r[[2]], quote=FALSE, row.names=TRUE);
}

And do the actual work:

meta_data <- read.table('./metadata.txt', sep = ' ', header = TRUE, stringsAsFactors = FALSE);
apply(meta_data, 1, process_row(r))
# You don't need a q() in an Rscript because it exits after the last line of code anyway.


Genomics Statistics Problem I have used a python script to identify target sequences in a DNA sequence file. There are two classes of sequence: coding and non-coding.

I have identified 728 sequences of interest. 597 of these fall into the coding regions and 131 of these fall into the non-coding regions. This is the equivalent of 18% non-coding. The total non-coding region in the sequence file is 13%.

Is there a statistical tool to demonstrate the python script identified target sequences in a non-random fashion?
If the script identified sequences that were randomly distributed then 13% of them would have been found in the non-coding region, from a total of 728 this seems like it should be reliable.

A Poisson or binomial would suffice.

I'm not at all sure the difference will significantly deviate from randomness.

Perl: 95% confidence interval of a non-parametric bootstrap for e.g. 1000 replicates is an alternative approach. I'd be more comfortable coding that solution in Perl than I would in Python, so is unlikely to be of much use here. Anyway, I would use a 2-dimensional array for your total data-set, the second dimension being coding/non-coding, I would then use the inbuilt random-number generator and multiply by the size of the array to randomly pluck out each element. 95% CI of the mean works by ranking the data and removing the top and bottom 2.5% and taking the interval for a single replicate (which is repeated 999 times more).

There will be a deep learning approach to this data, which would be easier to implement in Python, but I'd need to think about this because it is not simple classification.

Again my warning is that given the sample size, the result at best might be borderline. There are likely to be lots of interesting in your data, but coding/non-coding may not be one of them.


amino acids substitution probability I'm looking for the probabilities of amino acid substitutions in protein sequences.
A matrix of probabilities for every combination of amino acids (something like BLOSUM or PAM matrices but more universal and with raw probability values, not scores) is sought but alternatively, any method to ascertain this will be appreciated.
[EXTRA] If all diagonal values in such a matrix would be equal to 1, or you know any way to achieve that I'd be blessed.


How do I rename fasta headers? I am trying to use the seqkit replace command to replace chromosome names in the format chr_I, chr_II, ... to I, II, .... I am using the following command:
seqkit replace -p "(.)" --replacement "{kv}" --kv-file keyValues.txt mySequence.fasta

My keyValues.txt file contains the following:
chr_I   I
chr_II  II
...

The two columns are separated by a tabulation.
I get the following output:
[INFO] read key-value file: keyValues.txt
[INFO] 6 pairs of key-value loaded
[ERRO] pattern "(.)" matches multiple targets in "chr_I", this will cause chaos

I don't know what I am doing wrong. I chose the (.) pattern to match the whole header but it seems to be wrong. Any help would be appreciated.
Update
All the headers are shown below:
grep '>' mySequence.fasta
>chr_I
>chr_II
>chr_III
>chr_IV
>chr_V
>chr_X

Why not use sed instead?

sed -e 's/chr_I/I/' -e 's/chr_V/V/' -e 's/chr_X/X/' mySequence.fasta > mySeq.fasta

Or even simpler:

sed 's/chr_//' mySequence.fasta > mySeq.fasta


BAM to gene expression matrix (UMI counts per gene per cell),10X I am trying to reproduce some results of a scRNASeq experiment. However I am new to the server-side aspect of such analyses and am very confused at the moment.

The data provided by the authors of the paper is in .BAM format and from there I wish to derive a gene expression matrix (UMI counts per gene per cell). The authors stated that they did this using the 10X genomics pipeline, however it seems that .BAM outputs are sort of the "final version" of the pipeline and, the pipeline only really deals with truly raw sequencing files (bcl/fastq).

I have considered converting it back to fastq format and following the pipeline from the beginning, but it just seems like I would be moving backwards by doing such.

I have also found some R packages that can integrate BAM files in the environment like Rsamtools, but I would much rather do this on the server and then upload a csv into R for the downstream analysis as the files are +60GB each.

Essentially I am just asking for some advice, or a vignette/article I can read that would clear up how to proceed with the BAM outputs to derive the gene expression data. I have already read the cellranger pages on the 10X genomics website but they were of no help in answering this quesion.

Thank you for your time,

You should be able to parse out what you need using the tags in the .bam.  10xGenomics' website says what tags they add.

https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/output/bam

Going backwards would also involve parsing the tags, because you have to make two fastq files, and a simple bam -> fastq pipeline won't do that correctly.


PAML branch-site model: Identical lnL value between null and alternative models? I have been running PAML, via codeml, to try to identify positive selection in certain branches of my tree for one gene (but will eventually be trying to automate the analysis over multiple genes). This was done by setting model = 2 and NSsites = 2 in the codeml.ctl file. I attempted to run both the null and alternative models by setting fix_omega = 1 in the null.ctl and fix_omega = 0 in the alternative.ctl.

However, after running both models, they have returned exactly identical log-likelihood (lnL) values. I wasn't expecting to find anything significant, but the fact that they were identical struck me as odd. Out of curiosity, I altered the branch of interest in the newick tree file and reattempted the analysis, but again, the lnL values remain the same between the null and alternative models (although slightly different to the value returned with the different branch selected).

I have since repeated this many times, choosing a different branch(es) each time to see if I can get a difference between the null and alternative, but no luck. Am I missing something obvious? Should I really be getting the same exact result between null and alternative models, regardless of where in the tree I place #1? To reiterate, I wasn't expecting to find significant results, but was at least expecting (non-significantly) different numbers for lnL...

The issue with PAML is that it doesn't throw many exceptions (programming jargon for trapping and reporting errors) if something is wrong, generally the "exceptions" are identified in the failure to report an expected value in the output.

Possible explanation fix_omega is not turning the omega parameter on or off it is simply defining what the value of omega will be. If the default value, which is set at .4, but can be changed, is very close to the estimated value the likelihoods will be similar - or even identical.

Goto debugging As computer/data scientists, the very first thing we should do is just check there is run time variation between the options you are assessing via,

time codeml ./file

The thing to remember is fix_omega=0 is where Omega is being estimated from the data, whilst =1 means it will be supplied with a value. There should be longer run time on the "0" option.

Nitty gritty debugging

1 Please check you have

seqtype = 3

If not, then that is your answer, omega was never being set.

2 For debugging only make sure you have,

model = 1

This is global dN/dS for all the data not for a specific clade within the tree, the reason for this is given in the final step below

3, Then try a different NSsites model,

NSsites = 8

This should produce real differences. The NSsites options are,

  0:one w;1:neutral;2:selection; 3:discrete;4:freqs; * 5:gamma;6:2gamma;7:beta;8:beta&w;9:betaγ
  10:beta&gamma+1; 11:beta&normal>1; 12:0&2normal>1; * 13:3normal>0

w here is for omega

4 The proof of your system is that when you

fix_omega = 1
omega = .4 * try lots of different values, e.g. 0.3, 0.4, 0,7, 1, 1.5, 2

This should produce likelihood variation

Finally

model = 2

Then go to your tree file and make sure that the nomenclature used follows the manual, I think from memory the # is deployed. The output file should change to specify that you have two values for dN/dS from your treefile, thus the output from model=2 should vary from the output of model=1 if not the delineation has been ignored. If you do see two rates, great and now you are into parameter exploration. Parameter exploration is basically the options 1 to 4 described above, plus the stuff you mentioned.

Oh you mentioned using Newick format. I haven't used PAML for a while (dNdS stuff is ancient) but it it only used phylip format when I used it.  I suspect if it got the wrong tree format PAML will throw an exception, but it's worth checking.


Detecting differentially expressed genes with foldchange >= 2 and FDR < 0.05 I'm using edgeR for differential analysis. Using glmTreat function I'm detecting differentially expressed genes between Tumor and Normal. I have set the arguments like below:

tr <- glmTreat(fit, contrast=contrast.matrix, lfc=2)
tab2 <- topTags(tr,n=Inf, adjust.method = "BH")
keep <- tab2$table$FDR <= 0.05

This gave me differentially expressed genes with logFC > 2 at FDR < 0.05 and genes with logFC < 2 at FDR < 0.05. [There are few genes with logFC 1.9, 1.89 also]

I have also tried with lfc=log2(2).

tr <- glmTreat(fit, contrast=contrast.matrix, lfc=log2(2))
tab2 <- topTags(tr,n=Inf, adjust.method = "BH")
keep <- tab2$table$FDR <= 0.05

To my surprise I see DE genes with logFC >=1, FDR < 0.05 and also logFC <=1 at FDR < 0.05.

Here is an example of DE genes I also found when I used lfc = log2(2)

GeneSymbols logFC   unshrunk.logFC  logCPM  PValue  FDR
DPYD-IT1    1.14493675  1.560389373 0.67468395  0.002972841 0.018704644
LINC00945   1.144288473 1.525191847 0.693845996 0.00651175  0.038064312
ITPKB-AS1   1.14046752  1.484991704 0.703887505 0.006620218 0.038585862
AC105206.1  1.122495266 1.4871905   0.699167758 0.006545299 0.038194633
LINC02066   1.120660287 1.558458881 0.66016102  0.004344743 0.026424298
AC090985.1  1.10235045  1.58968763  0.64343419  0.002745089 0.017426607
KCNJ6-AS1   1.100399688 1.717028637 0.613324799 0.002224039 0.01438451
LINC01494   1.090270138 1.490899364 0.676937984 0.006770608 0.039348027
AC099541.1  1.080957092 1.53817727  0.651356989 0.008449298 0.048356352
AL133396.2  1.012984207 1.512530231 0.63130643  0.008090259 0.046411928
AC011363.1  1.01207753  1.499708945 0.635740622 0.007226791 0.041797256

I'm very confused with these foldchange cutoffs.

Which one should I use to select differentially expressed genes based on fold change >=2 and FDR < 0.05?

The edgeR authors recommend that you use a relatively low logFC threshold for glmTreat such as lfc=log2(1.2). A lfc value as high as lfc=log2(2) is seldom required and only would only be appropriate for datasets with really large amounts of DE.
The first threshold you used of lfc=2 is equivalent to lfc=log2(4) and is too high for any dataset.
There's a few subtleties here:

First, there's the difference between logFC and FC as explained by Devon.

The lfc in glmTreat is not a fold-change cutoff! It is instead a testing threshold. All the significant genes will have actual logFCs well above the threshold.

The lfc threshold applies to the raw fold-changes (unshrunk.logFC in the output), not to the shrunk versions (logFC in the table).

The intent of glmTreat is to give a better gene ranking than would be provided by FDR or fold-change, so please don't re-sort the results by fold-change -- that would tend to defeat the purpose of the method.

I understand the documentation for glmTreat is a bit light, but you can get an idea of how to use it from our workflows.

https://bioconductor.org/packages/release/workflows/vignettes/RnaSeqGeneEdgeRQL/inst/doc/edgeRQL.html

https://bioconductor.org/packages/release/workflows/vignettes/RNAseq123/inst/doc/limmaWorkflow.html

I understand the simple fold-change cutoff of FC > 2 with FDR < 0.05 is often used in the literature, but the edgeR authors do not recommend it and hence we don't provide an automatic way to get it from the edgeR pipeline.
For RNA-seq, fold-change cutoffs tend to prioritize genes with low expression and low counts.
The glmTreat approach is more subtle but we think gives better results, both in terms of biological relevance and statistical correctness.


How to convert txt.gz into fastq.gz? I'm a newbie to bioinfo. I'm trying to use galaxy.networkanalyst.ca for the analysis of RNA-Seq data and for that I need to convert txt.gz files into fastq.gz files.

Link: GSM4273445 is one sample and the whole series is from here.
I've Googled to find out if there are online tools to convert it, but there aren't any. Besides, there is no discussion addressing this format conversion on researchgate.
I intend to use this RNA-Seq data and investigate the relative expression of some specific transcripts of interest. I have been introduced to 'galaxy.networkanalyst.ca' where the very first step is to upload fastq files. In addition to this I was aware of ncbi GEOdatasets as source of open RnA-Seq data.
Can anyone suggest a method or tool for this?

I've had a quick look at the file you linked to and this is a featureCounts file - this list the number of reads for each gene/transcript (COUTNS.txt in the file name pretty much gave it away that it was going to be a file of, well counts).
I can't help you with Galaxy, but if you want to do an RNA-seq analysis, the count data + metadata is likely all you'll need.
If you really want to do the preprocessing and mapping yourself, the fastq files can be downloaded from the SRA archive (SRP242087) - click on the SRA link in your whole series link.


Alignment QC differences between HISAT2 and Qualimap I used hisat2 for aligning reads to to the genome. I have an alignment summary for sample1 as follows:

89038751 reads; of these:
  89038751 (100.00%) were paired; of these:
    5641019 (6.34%) aligned concordantly 0 times
    67665552 (76.00%) aligned concordantly exactly 1 time
    15732180 (17.67%) aligned concordantly >1 times
    ----
    5641019 pairs aligned concordantly 0 times; of these:
      651415 (11.55%) aligned discordantly 1 time
    ----
    4989604 pairs aligned 0 times concordantly or discordantly; of these:
      9979208 mates make up the pairs; of these:
        5887822 (59.00%) aligned 0 times
        2656902 (26.62%) aligned exactly 1 time
        1434484 (14.37%) aligned >1 times
96.69% overall alignment rate

To check the alignment quality I used Qualimap. In the report I see the following:

Reads alignment

Number of mapped reads: 172,189,680
Total number of alignments: 244,302,519
Number of secondary alignments: 72,112,839
Number of non-unique alignments:    22,973,098
Aligned to genes:   2,081
Ambiguous alignments:   173
No feature assigned:    178,975
Not aligned:    5,887,822

Reads genomic origin

Exonic: 2,081 / 1.15%
Intronic:   6,866 / 3.79%
Intergenic: 172,109 / 95.06%
Intronic/intergenic overlapping exon:   302 / 0.17%

Transcript coverage profile

5' bias:    0.69
3' bias:    2.85
5'-3' bias: 0.39

Can anyone tell me why the numbers look different? Did I go wrong somewhere? Help me in understanding the results.

89038751 pairs times 2 is 178 million reads. That times the approximately 97% alignment rate is around 172 million reads, which is what qualimap is reporting.

A "non-unique" alignment is one that can map to multiple places (either equally well, or where all of the alignments are of similar enough quality). Such alignments have a single "primary" alignment and one or more "secondary" alignments, which is where the "secondary" metrics come from.

For the qualimap output, there are additional metrics describing how reads overlap genes (or not). Of note are the "ambiguous" alignments, which are those overlapping multiple genes. "No feature" metrics, then, refer to alignments not overlapping any gene (see the qualimap documentation for further details).

My guess is that you have such a high "intergenic" rate because you've either used the wrong GTF file or have told qualimap the wrong strand setting. For unusual organisms one expects many reads to arise from unannotated genes. This is not the case for common model organisms, such as mice and humans, where there are very few remaining commonly expressed unannotated genes.


How to subset samples from a VCF file? I have VCF files (SNPs & indels) for WGS on 100 samples, but I want to only use a specific subset of 10 of the samples. Is there a relatively easy way to pull out only the 10 samples, while still keeping all of the data for the entire genome?

I have a script that allows me to pull out regions of the whole genome for all 100 samples, so if I could do something similar but only put regions for the 10 samples that I want that would be ideal.

Thanks for the help!

Bcftools has sample/individual filtering as an option for most of the commands. You can subset individuals by using the -s or -S option:

  -s, --samples [^]LIST

  Comma-separated list of samples to include or exclude if prefixed with "^". Note that in general tags such as INFO/AC, INFO/AN, etc are not updated to correspond to the subset samples. bcftools view is the exception where some
  tags will be updated (unless the -I, --no-update option is used; see bcftools view documentation). To use updated tags for the subset in another command one can pipe from view into that command. For example:

  -S, --samples-file FILE

  File of sample names to include or exclude if prefixed with "^". One sample per line. See also the note above for the -s, --samples option. The command bcftools call accepts an optional second column indicating ploidy (0, 1 or 2) or sex (as defined by --ploidy, for example "F" or "M"), and can parse also PED files. If the second column is not present, the sex "F" is assumed. With bcftools call -C trio, PED file is expected. File formats examples:

sample1    1
sample2    2
sample3    2

  or

sample1    M
sample2    F
sample3    F

  or a .ped file (here is shown a minimum working example, the first column is ignored and the last indicates sex: 1=male, 2=female):

ignored daughterA fatherA motherA 2
ignored sonB fatherB motherB 1

Example usage:

bcftools view -s sample1,sample2 file.vcf > filtered.vcf
bcftools view -S sample_file.txt file.vcf > filtered.vcf

See the bcftools manpage for more information.


How to combine two Genome-wide Association Study (GWAS)? I did a GWAS analysis in the past for antibiotic resistance of E. Coli and the results were interesting (matching the literature). I did a new GWAS analysis for some new samples, but the results are not in agreement with the previous GWAS and the literature. So, I was wondering if the combination of them would be in agreement with the literature, as having more samples to support an argument is always a good idea.

So, my question is how I can combine two GWAS study and what quality control steps do I need to follow?

The first way is to combine two samples and re-run GWAS. A second way is to do a meta-analysis using GWAS results.


Standalone Blast+; Automating searches, Blastn formatting, interpretation I am struggling with using NCBI's standalone blast, particularly the blastn feature.

I have a large number of nucleotide sequences that I want to use blastn to ID. The goal is to retrieve information from the blast program and eventually import this information to excel. I want to know the capabilities of, or the easiest way, to format the output. I am running a remote search to NCBI's nucleotide database and want to filter and organize the information retrieved as much as possible. The ultimate goal is automate a way to run a remote search to NCBI's nucleotide database, and to work with the results in excel.

I am extremely new to this and struggling to find the best path. From what I can gather, many prefer using external programs, scripts, parsers, interpreters etc to organize their results. I am okay with doing this as well, but would still like to format the output of my blastn search as much as possible beforehand. It is returning a lot of data to me that I don't need.

Can I automate this in excel, and just use VBA/macros to further organize my data? It seems some people prefer to just import the output text file to a spreadsheet and work from there. I have also seen options like BioEdit, eutilities, APIs... I'm not sure which route I should invest my time into figuring out.

Thanks so much for lending your time and expertise!

You just have to change the output format using the -outfmt parameter. To have an output in a tabular format just use -outfmt 6. The resulting table have 12 columns containing the following informations queryacc subjectacc pident length mismatch gapopen querystart queryend subjectstart subjectend evalue bitscore. You may restrict the number of hits per query using -max_target_seqs 5 for example.


Get results of keyword search on Pfam via python script I'm interested in all proteins that are in any way associated with Danio rerio. I decided to look them up at Pfam data base and when I just make a keyword search, I get a a nice list which looks like this http://pfam.xfam.org/search/keyword?query=danio+rerio&submit=Submit.
Important thing to note is that most of the associations of these protein families can be found in seq_info column.

My question is: how to extract that list from a python script? Moreover, I want to make a script that when run, will fetch the latest release, filter out the proteins that cannot be found in D. rerio and do some additional processing.

I've been searching for ways to directly download files from Pfam ftp server and mine for the proteins that are associated with D. rerio, however, I didn't find a single way to do it. I simply cannot see where do they store the information about species in which the protein can be found.

I would like to retrieve everything about protein families of interest, because the next step I will need to do is to find the DNA binding proteins and for finding a way to filter out the ones of no interest, I would liketo have as much data as possible.

If I understand your question correctly, then you can use the Pfam proteom search to get the data you are looking for.

Go to the

  Pfam website -> Browse -> select the letter D from proteoms category
  -> D. rerio

There you can find a download link under the section Domain composition:

  You can download a list of all regions for this proteome from our FTP
  site.

The download link seems to be bound to the species ID and if you want to automate your data retrieval you may be able to simply download the file from there in your python script.


How to count RNA sequencereads using custom made python scripts? I am trying to do RNA seq analysis and my goal is to filter gene counts less than 5 using custom made python scripts. The code chunk goes as follows
#!/usr/bin/python
#Load relevant modules
import sys

#Read in the HTseq gene count output files for each of the 8 samples
inputlist = [sys.argv[1], sys.argv[2], sys.argv[3], sys.argv[4], sys.argv[5], sys.argv[6], sys.argv[7], sys.argv[8]]
##Determine the genes with sufficient coverage
#Set up a variable to hold genes with sufficient coverage in any of the 8 samples
genelist = []

#For each HTseq file
for i in range(len(inputlist)):
    #Open it
    inputfile = open(inputlist[i], 'r')
    #For each line
    for line in inputfile:
        #Determine the number of reads/transcripts for this gene
        splitline = line.strip().split('\t')
        #Does it pass the threshold?
        if int(splitline[1]) >= 5:
            #If so, add it to list
            if not splitline[0] in genelist:
                genelist.append(splitline[0])

    #Close the file
    inputfile.close()

#Write out the list of sufficiently expressed genes
outputlist = open(sys.argv[9], 'w')

for gene in genelist:
    outputlist.write(gene + '\n')

outputlist.close()
##Filter each of the HTseq files for only sufficiently expressed genes across our samples
#For each HTseq file

for i in range(len(inputlist)):
    #Open it
    inputfile = open(inputlist[i], 'r')

#Open an output 'filtered' file
outputfile = open(sys.argv[i+10], 'w')

#For each line
for line in inputfile:
    #Determine the gene
    splitline = line.strip().split('\t')

    #Is it in our list?
    if splitline[0] in genelist:
        #If so, write out the original line
        outputfile.write(line)

#Close the files
inputfile.close()
outputfile.close()

python filtercounts.py 'wth1count.txt' 'wth2count.txt' 'wth3count.txt' 'wth4count.txt' '99h1count.txt' '99h2count.txt' '99h3count.txt' '99h4count.txt' 'genelist.txt' 'wth1filtered.txt' 'wth2filtered.txt' 'wth3filtered.txt' 'wth4filtered.txt'  '99h1filtered.txt' '99h2filtered.txt' '99h3filtered.txt' '99h4filtered.txt'
The program works fine but at the end of the code execution I am not getting the last two output files '99h3filtered.txt' '99h4filtered.txt'.

What software are you using for the next step?  DESeq in R will do this filtering for you in one line.
I think a more typical way to do this would be to make your entire gene count matrix, and then remove rows where the max value for that row is < 5.  Your way, you have to go through each file twice; once to know which genes to keep, and then again to actually remove lines.


Bioinformatics approach to dentifing potential PCR primer sequences for transcribed gene I have an annotated transcriptome and would like to develop PCR primers for particular transcribed genes. My species is a non-model plant. Can I use BLAST or another tool to identify potential PCR primer sequences?  Or more generally, is there a workable bioinformatics approach to identifying potential PCR primers for transcripts arising from a non-model organism? Any literature refs would be great too.

It depends if you want to amplify your target from genomic DNA or cDNA.

Depending on intron size, primers that anneal to different exons may not be able to amplify across the intron(s), and if you're unlucky enough to design a primer that spans a splice junction, you'll definitely be stuffed if you're working with genomic DNA. I would try blastn or translate to an amino acid sequence and use tblastn to see if you can identify the intron/exon structure in a related species. If you don't have a genome sequence available for anything close, I think your only option is to target as small an amplicon as possible and hope for the best.

If you are working with cDNA, any standard primer design tool should do


How to calculate the number of cells in different cell clusters with Seurat v3? I was supposed to use the code below to get the number of cells per cluster. The code works with Seurat version 2, but while using version 3 I got the error no slot of name "ident" for this object of class "Seurat". Does anyone have any idea about how to deal with this?

# Calculate number of cells per cluster from object@ident

cell.num <- table(shox2crepos_mtmg@ident)

# Add cell number per cluster to cluster labels

ClusterLabels = paste("Cluster", names(cell.num), paste0("(n = ", cell.num, ")"))

# Order legend labels in plot in the same order as 'ClusterLabels'

ClusterBreaks = names(cell.num)

# Plot tSNE with new legend labels for clusters

TSNEPlot(object = shox2crepos_mtmg, do.return = T) +
  scale_colour_discrete(breaks = ClusterBreaks,
                        labels = ClusterLabels) +
  labs(x = "t-SNE 1",
       y = "t-SNE 2")

The "identities" can be accessed with the Idents() function in Seurat v3.

You should change cell.num <- table(shox2crepos_mtmg@ident) to

cell.num <- table(Idents(shox2crepos_mtmg))


What motif finding software is available for multiple sequences ~10Kb? I have around ~3,000 short sequences of approximately ~10Kb long. What are the best ways to find the motifs among all of these sequences? Is there a certain software/method recommended?

There are several ways to do this. My goal would be to:

(1) Check for motifs repeated within individual sequences

(2) Check for motifs shared among all sequences

(3) Check for the presence of "expected" or known motifs

With respect to #3, I'm also curious if I find e.g. trinucleotide sequences, how does one check the context around these regions?

Thank you for the recommendations/help!

The MEME Suite web site contains a collection of tools for motif analysis (I'm one of the maintainers). It contains two de novo motif discovery tools: MEME and DREME. Public web applications are provided, but you can also download and build command line tools for a local installation.

For your first goal you could use MEME and select the "Any number of repetitions model" (ANR). For your second goal, you'd use MEME with the "Zero or One Occurrences Per Sequence" (ZOOPS) model. For your third goal you could use FIMO (Find Individual Motif Occurrences), and one or more of the motif databases provided on the software and database download page.

It sounds like your sequence data is about 30Mb. The MEME web application is limited to 60kb of sequence data, so you'd have to install a local copy of the MEME Suite. MEME would take a long time to analyze a 30Mb sequence database unless you have MPI configured, and lots of cores available. You might want to consider analyzing a randomly selected subset of your sequences. The running time of MEME grows as the cube of the number of sequences.

For short motifs, you may want to use DREME rather than MEME. DREME is better than MEME at identifying short motifs, but is limited to motifs <= 8 positions wide.


gene-level versus transcript-level analysis Traditionally, RNA-seq data was quantified on gene level. Newer methods quantify on transcript/isoform level. For example, Kallisto only outputs transcript-level abundances. From the DESeq2 vignette:

  A newer and recommended pipeline is to use fast transcript abundance
  quantifiers upstream of DESeq2, and then to create gene-level count
  matrices for use with DESeq2 by importing the quantification data
  using the tximport package.

Is it just for consistency and/or simplicity that the values are converted to gene level? Would it better to proceed with transcript-level table or does that violate some DESeq2 (or similar tools) assumptions?

There are a variety of reasons people use gene-level quantitations.

Transcript-level differences are difficult to biologically interpret. Let's be honest, few groups are likely to put in the work required to determine what these might mean. Most genes have at least some level of characterization, so it's much more biologically tractable to think in terms of gene-level changes. Granted, "genes" are largely artificial constructs and in an ideal world one would always think in terms of transcripts, since they're what's really changing. But we don't live in an ideal world.
Statistical power increases with counts, so if you're splitting your counts among a handful of isoforms you could be missing something. In practice this isn't a big problem, since there are typically obvious major isoforms, but that won't always be the case.
"Mo transcripts mo problems", to channel the bioinformatics notorious B.I.G. It's common practice for people to use a subset of transcripts when using tools like Kallisto and Salmon, since the results usually end up being better. But which subset of transcripts should you use? Those in refseq? Gencode basic? A subset of one of those? Do we really understand all the consequences of choosing one over the others? (No, we don't) We all have our own preferences, but this isn't a solved problem.