BioMart Ensembl id imperfect retrieval I have a list of gene symbols:

c("cd45", "Tmem119", "CD11b", "P2Yr12", "Csf1r", "Bst2", "Cd74",
"Cx3cr1", "Trem2", "Lyz2", "GLAST", "GFAP", "ALDH1L1", "Aldoc", "Aqp4",
"Glul", "S100a", "Olig1", "Olig2", "Olig3", "Mbp", "Pdgfra", "Pecam",
"Cldn5", "Cldn10", "Epas1", "Crip1")

If I feed it to BioMart to get ensembl ids:

mart <- useDataset("mmusculus_gene_ensembl", useMart("ensembl"))
list <- getBM(filters= "mgi_symbol", attributes= c("ensembl_gene_id",
                           "mgi_symbol","description"),values=symbols,mart= mart)

the following genes are missed:

"cd45"   "cd11b"  "p2yr12" "glast"  "s100a"  "pecam"

All of them are pretty well-known genes, and I can manually find their ensembl ids by googling it, for instance:

http://www.informatics.jax.org/marker/MGI:97810

I tried supplying aliases, but it does not change the output. So, from my understanding BioMart is not working properly either because I am doing something wrong here or because BioMart itself is not a good tool to use. Is there a better way of getting the mapping that would map all of the gene symbols?

The genes that are missed are probably not official mgi symbols.

You might wanna look them up at mgi:

cd45 -> Ptprc

cd11b -> Itgam

p2ry12 -> P2ry12 (?) (CAPS sensitive?)

glast -> Slc1a3

s100a -> S100a1

pecam -> Pecam1 (?)

I suspect biomaRt is not the problem here, the only better way I can come up with is to download the ENSEMBL annotation (GTF) file, and get the symbol out of this annotation with AWK or in R.

Here an example in R, since you don't want to use AWK:

file <- read.table("Mus_musculus.GRCm38.78.gtf", sep = "\t", comment.char = "#")

exons <- subset(file, V3 == 'exon')

V9 <- file[,"V9"]

ensembl_names <- sub("gene_id ","",V9)

split <- strsplit(ensembl_names, ";")

splitting <- sapply(split, function(x) x[1])

split_symbol <- strsplit(ensembl_names, "gene_name ")

splitting_symbol <- sapply(split_symbol, function(x) x[2])

splitting_symbol2 <- strsplit(splitting_symbol, "; ")

symbol_names <- sapply(splitting_symbol2, function(x) x[1])

ensemble_symbol_big <- cbind(splitting, symbol_names)

ensembl_symbol <- unique(ensemble_symbol_big)

colnames(ensembl_symbol) <- c("ENSEMBL", "Symbol")

write.table(ensembl_symbol, "ensembl_symbol.txt", sep = "\t", row.names = F)


Hosting IGV web on a server I have my analysis on a cloud server. The data is quite large, so copying it locally is not really an option. I want to be able to host IGV as a server on that server, so that it can access those files.

I installed IGV web (https://github.com/igvteam/igv-webapp) on my server, and ran it, but the Tracks → Local File option opens a file browser on my local machine, it doesn't let me select files on the server.

Is there any easy way to host IGV web on a cloud server?

You can use X11 forwarding to run IGV on the headless server where the data is located and simply forward the GUI to your client machine using the X11 protocol. This is a solution I use consistently for my research when working with 20+ tracks that I want to keep on the server.

Download and install IGV for linux from the IGV site on to your server.
Connect to the server using ssh with X11 forwarding enabled

ssh -X user@server
Navigate to the recently downloaded and unzipped IGV folder. Execute the igv.sh script.

./igv.sh

The script launches the IGV session with 4gb of memory. Edit the script to launch with more, 32gb in this edited excerpt.

exec java -showversion --module-path="${prefix}/lib" -Xmx32g \
@"${prefix}/igv.args" \
-Dapple.laf.useScreenMenuBar=true \
-Djava.net.preferIPv4Stack=true \
--module=org.igv/org.broad.igv.ui.Main "$@"


Error while Importing fastq sequences to jMHC software I am trying to analyze my .fastq files (MiSeq exon reads) with the jMHC software, but output files I have obtained have no results.

I think I might be doing something wrong during the first step: "connect to database".
Because on the "import/extracting" step, I open all my .fastq files, click "Start", and it shows this message:

  "Success! Processed: 139, Imported: 0".

So I assume it is reading the files OK, but not importing them to the "database"? Could it be the extension of my files? The software says it operates on FASTA files. Are .fastq files not compatible?

Anybody that has successfully run this program could give me a hand?

Yes, use .fasta , you can try to convert your fastq files with NextGen Workbench for example.


How to retrieve InChI key for KEGG compound? I would like to retrieve the InChI representation for a given KEGG compound but I fail to find a direct solution for this (preferably doing it via bioservices).

One could do it via ChEBI like this:

from bioservices import *

kegg_con = KEGG()

kegg_entry = kegg_con.parse(kegg_con.get('C00047'))

chebi_con = ChEBI()

chebi_entry = chebi_con.getCompleteEntity('CHEBI:' + kegg_entry['DBLINKS']['ChEBI'])

print chebi_entry.smiles
print chebi_entry.inchi
print chebi_entry.inchiKey

which will print

NCCCC[C@H](N)C(O)=O
InChI=1S/C6H14N2O2/c7-4-2-1-3-5(8)6(9)10/h5H,1-4,7-8H2,(H,9,10)/t5-/m0/s1
KDXKERNSBIXSRK-YFKPBYRVSA-N

Is there a direct way which does not require to go via ChEBI?

I'm not aware of a functionality in the KEGG service that would allow you to do that directly. I believe your solution (using e.g. ChEBI) is the correct one.

You may have other useful service such as unichem, chemspider (you would need a login though for that one) could help as well.

Disclaimer: I'm the bioservices main author


Are there certain alignment methods/tools which perform better with a high density of indels? I have a set of experiments which should result in WGS reads with a high density of indels.

Question: Are there certain alignment tools/methods which perform respectively "better" with these reads? I'm not very interested in performance, but rather "better" in terms of dealing with the presence of many deletions/insertions.

Options I know of:

BWA  (I think there are special settings whereby I could run bwa
for this case, e.g. see this)
BBMap

Minimap2 attempts to solve a lot of problems that are not properly dealt with by other mapping software. It's written by the same author as BWA, but designed to deal with situations where the reference and query sequences are quite different (e.g. for noisy long reads).


Why sequence the human genome at 30x coverage? A bit of a historical question on a number, 30 times coverage, that's become so familiar in the field: why do we sequence the human genome at 30x coverage?

My question has two parts:

Who came up with the 30x value and why?
Does the value need to be updated to reflect today's state-of-the-art?

In summary, if the 30x value is a number that was based on the old Solexa GAIIx 2x35bp reads and error rates, and the current standard Illumina sequencing is 2x150bp, does the 30x value need updating?

The earliest mention of the 30x paradigm I could find is in the original Illumina whole-genome sequencing paper: Bentley, 2008. Specifically, in Figure 5, they show that most SNPs have been found, and that there are few uncovered/uncalled bases by the time you reach 30x:

These days, 30x is still a common standard, but large-scale germline sequencing projects are often pushing down closer to 25x and finding it adequate. Every group doing this seriously has done power calculations based on specifics of their machines and prep (things like error rates and read lengths matter!).

Cancer genomics is going in the other direction. When you have to contend with purity, ploidy, and subclonal populations, much more coverage than 30x is needed. Our group showed in this 2015 paper that even 300x whole-genome coverage of a tumor was likely missing real rare variants in a tumor.

On the whole, the sequence coverage you need really depends on what questions you're asking, and I'd recommend that anyone designing a sequencing experiment consult with both a sequencing expert and a statistician beforehand (and it's even better if those are the same person!)


Where to download JASPAR TFBS motif bed file? I am interested in determining if any transcription factor binding site motifs are enriched in some BED files from a DNA methylation experiment.

I am looking for a database that has BED Files containing regions enriched for specific transcription factor binding motifs in homo sapiens.

I know that JASPAR contains TFBS motifs, but I can't find a download for regions of hg19 that are enriched with those motifs.

EDIT:  I suppose I could determine TF motif enriched regions using HOMER if there are no curated databases with this information, but I would rather not.  Your help is appreciated!

You could make your own such file with FIMO, the JASPAR MEME files available via the MEME download site, and your genome build of interest (per-chromosome FASTA files from UCSC goldenpath, say), e.g., for assembly hg19 and UCSC chromosome naming scheme:
$ for chr in `seq 1 22` X Y; do echo ${chr}; wget -qO- http://hgdownload.cse.ucsc.edu/goldenpath/hg19/chromosomes/chr${chr}.fa.gz | gunzip -c - > hg19.chr${chr}.fa; done

For each chromosome (omitting the for loop code for brevity), you can use the fasta-get-markov in the MEME suite to generate a 5th-order Markov model of kmer frequencies:
$ fasta-get-markov -m 5 hg19.chr${chr}.fa > hg19.chr${chr}.5th_order_background_model.txt

Here is one way to download a non-redundant, vertebrate JASPAR database (MEME-formatted):
$ wget http://jaspar.genereg.net/download/CORE/JASPAR2020_CORE_vertebrates_non-redundant_pfms_meme.txt

To run a FIMO query with the 5th-order background model and JASPAR database, using a p-value threshold of 1e-4:
$ fimo --verbosity 1 --bgfile hg19.chr${chr}.5th_order_background_model.txt --thresh 1e-4 --text JASPAR2018_CORE_vertebrates_non-redundant_pfms_meme.txt hg19.chr${chr}.fa \
  | tail -n+2 \
  | awk -vOFS="\t" -vFS="\t" '{ print $3, $4, $5, $1":"$2">"$10, $8, $6 }' \
  | sort-bed - \
  | starch --omit-signature - \
  > hg19.chr${chr}.bed.starch

The output of a FIMO search is a text file which can be modified into a sorted BED file via sort-bed, an example of which is shown above.
To get a faster answer, each per-chromosome query can be put onto its own node on a computational cluster.
Using a p-value threshold of 1e-4 will result in very large files (when uncompressed). Using tools to compress BED files is useful here. Use of BEDOPS starch to compress BED files (as shown above, for instance) will give better compression performance than bgzip with per-chromosome random access and integration with BEDOPS bedops and bedmap CLI tools. The bgzip option offers integration with tabix and offers finer-grained random access, but it will use more disk space and doesn't integrate with BEDOPS tools.


estimate genome size: kmer-based approach from PacBio reads Can anyone suggest a software/method for kmer analysis using PacBio reads (RSII)?

Something similar to Jellyfish, that I saw in a nice tutorial - but must be suitable for long, noisy reads. kmercounterexact from BBMapTools might also be an option, but again tricky with long, noisy reads.

The closest solution I could find was to estimate the genome size from overlap statistics (for the raw reads, the error corrected preassembled reads or by mapping the raw reads to the assembled contigs) by:

running an initial assembly using a small seed read length
then plot the pre-assembled read overlap histogram.

What do you think of this approach?

I would like to use this kmer analysis to estimate genome size, but for now I could not find much. I do not have companion Illumina reads - if so, I could try the approach proposed by Carvalho et al. 2016

I have reads & in-house genome assemblies (Canu) from different isolates of the same species of fungus to compare.

I don't think there is a method that would estimate a genome size using raw long reads.

The genome size estimates based on raw reads are done by fitting a model to kmer spectra (for instance Genomescope). The kmer spectra built from long reads are really messy due to the high error rate of long reads. That makes fitting of a model quite difficult. These methods assume small error rates. If the error rate is as high as 0.1 (which is normal for long reads), then the probability of having one kmer of 23 bases correct is ~8%. Meaning that if the coverage is 100x, just ~8x of homozygous kmers would be the real genomic kmers and it's very likely that kmer peaks of the sequencing errors, heterozygous kmers and homozygous kmers will just blend into one.

I can imagine, that if your coverage would crazy high (like 600x or more) the kmer spectra could separate errors and genomic kmers. You could also try to decrease kmer size to increase the proportion of correct kmers.

Another option would be to use a trick that they used to generate solid kmers for genome assembly using generalized de brujin graph using noisy reads. Maybe it would be better to take a look on the newer version of the assembler, that is called Flye.

But to be honest, I think for a long reads, a genome assembly is much more reliable way how to predict a genome size than any of the kmer tricks. I would simply assemble the genome, remap back reads, check the uniformity of coverage and conclude the genome size out it.

-- edit --

I tried to find some comprehensive evaluation of the assembly sizes compared to genome sizes, but I could not, so I just checked two examples, one big, one small and they both match quite well. Ammopiptanthus nanus was assembled to 823.74 Mb while kmer estimate was 889 Mb (ref). Plasmodium falciparum was assembled to 23.5 Mb while the original genome size estimate is 22.8 Mb. I know it's just two, but there is also a good reason to think that the genome size will be reasonable, feel free to add here more examples.

If you would be worried about contamination in the assembly, you can first run it though one of the methods to detect a contamination (like Blobology).


How to compare transcriptomic profiles of two cell types (single cell RNA-seq)? I found this interesting Single RNA-seq data set in GEO, but I am not sure how to analyze it appropriately.

They have deposited transcriptomic profiles of human and mouse pancreatic islets (pancreatic cells: Beta cells, Delta, etc).  The problem I see is that the different pancreatic cell types are not in equal numbers: 872 Beta cells and 214 Delta cells.

I want to see whether the expression of a set of genes are correlated between the two cell types. I don't think that I can visualize the correlation between the expression of the Beta and Delta cells on a scatter plot if they are unequal in number. Is there a workaround to include all the cells in the analysis, without losing good quality transcriptome of beta/delta cell?

What I would like to do is to compare the expression of two pancreatic cell types (Beta v.s. Delta cells) using a scatter plot. I have a set of genes that I am interested in, and would like to see their expression profile for beta and delta cells. As far as I know, it is only possible to visualize this using a scatterplot (X:Beta vs Y:Delta) where the number of isolated beta cells is the same as the number of delta cells.

Question: Given unequal number of isolated panacriatic cells, what would be an appropriate way to compare the expression profiles? Should I just ignore the extra ones?

Any idea how to approch this problem?

Thanks,

You are looking for a way to compare expression profiles between cell types, and your data is counts of genes expressed per cell. Your issue is that the data is highly dimensional. It has as many dimensions as you have genes. To be able find useful information, you must reduce these dimensions using a statistical technique such as Principal COmponents Analysis (PCA), t-Distributed Stochastic Neighbor Embedding (t-SNE), or Multidimensional Scaling (MDS). The result of these analyses is a reduction of your thousands of dimensions to vectors (eigenvectors in PCA) of variance which you can plot against each other in a scatter plot.

I have been working with Single Cell RNA-seq data and have been exploring the difference in gene expression in a certain type of cells between timepoints. Seurat (R package) is a very powerful tool for scRNA-Seq data analysis.

Here is a a tSNE of my data. I have added the timepoint labels, but you can see that after performing the dimensionality reduction the two different timepoints cluster by their gene expression profiles. I can then use Seurat to see which specific genes are differentially expressed between the clusters and a whole bunch of other fun things.


Paired end sequencing with R1 and R2 of different length. Possible? There are commercial sequencing kits/sequencers that allow for paired end sequencing in which the two reads obtained for each fragment are of different length?

For example 150bp for R1 and 50bp for R2?

Is it at least possible in principle? There are some published speculations about that possibility?

As indicated in the comments, yes, you can definitely tell standard Illumina sequencers to sequence mates in a pair to different lengths. This is quite common in single-cell RNA-seq, where R1 is often just a cell barcode and an UMI, so sequencing it to 150 bases would be a huge waste or reagent. One caveat to this is that you can't buy reagents sized for these custom length, just things like 150 base PE or other standards. So unless you end up saving and pooling reagents after such runs you end up paying the same for projects with one shorter read than if you'd just directly sequenced everything to the same length. Given that sequencer runs only take a day or two at the moment, this also doesn't result in that immense of a time savings.


Tool or script to parse annotated VCF files I've annotated VCF files using snpEff and looking for a tool or script to parse the VCF file and clean up the file to make it interpretable for a biologist.


cooler out of memory I ran into the following error when using cooler to load a contact pairs file:

ubuntu@ip-172-31-8-42:/data/hg-tmp$ FILE=GSE109344_monkey_fibro_allValidPairs.txt.gz; sudo docker run -v $(pwd):/data $TAG cooler cload pairs -c1 2 -p1 3 -c2 5 -p2 6 /data/chromInfo.rheMac2.
simple.txt:${BINSIZE} /data/$FILE /data/${FILE}.cool;
Traceback (most recent call last):
  File "/usr/local/bin/cooler", line 11, in <module>
    sys.exit(cli())
  File "/usr/local/lib/python2.7/site-packages/click/core.py", line 764, in __call__
    return self.main(*args, **kwargs)
  File "/usr/local/lib/python2.7/site-packages/click/core.py", line 717, in main
    rv = self.invoke(ctx)
  File "/usr/local/lib/python2.7/site-packages/click/core.py", line 1137, in invoke
    return _process_result(sub_ctx.command.invoke(sub_ctx))
  File "/usr/local/lib/python2.7/site-packages/click/core.py", line 1137, in invoke
    return _process_result(sub_ctx.command.invoke(sub_ctx))
  File "/usr/local/lib/python2.7/site-packages/click/core.py", line 956, in invoke
    return ctx.invoke(self.callback, **ctx.params)
  File "/usr/local/lib/python2.7/site-packages/click/core.py", line 555, in invoke
    return callback(*args, **kwargs)
  File "/usr/local/lib/python2.7/site-packages/cooler/cli/cload.py", line 477, in pairs
    map(pipeline, reader),
  File "cytoolz/functoolz.pyx", line 493, in cytoolz.functoolz.Compose.__call__
  File "/usr/local/lib/python2.7/site-packages/cooler/create/_ingest.py", line 425, in _aggregate_records
    .aggregate(agg)
  File "/usr/local/lib/python2.7/site-packages/pandas/core/groupby/generic.py", line 1315, in aggregate
    return super(DataFrameGroupBy, self).aggregate(arg, *args, **kwargs)
  File "/usr/local/lib/python2.7/site-packages/pandas/core/groupby/generic.py", line 186, in aggregate
    result, how = self._aggregate(arg, _level=_level, *args, **kwargs)
  File "/usr/local/lib/python2.7/site-packages/pandas/core/base.py", line 498, in _aggregate
    result = _agg(arg, _agg_1dim)
  File "/usr/local/lib/python2.7/site-packages/pandas/core/base.py", line 449, in _agg
    result[fname] = func(fname, agg_how)
  File "/usr/local/lib/python2.7/site-packages/pandas/core/base.py", line 432, in _agg_1dim
    return colg.aggregate(how, _level=(_level or 0) + 1)
  File "/usr/local/lib/python2.7/site-packages/pandas/core/groupby/generic.py", line 760, in aggregate
    return getattr(self, func_or_funcs)(*args, **kwargs)
  File "/usr/local/lib/python2.7/site-packages/pandas/core/groupby/groupby.py", line 1227, in size
    result = self.grouper.size()
  File "/usr/local/lib/python2.7/site-packages/pandas/core/groupby/ops.py", line 225, in size
    ids, _, ngroup = self.group_info
  File "pandas/_libs/properties.pyx", line 36, in pandas._libs.properties.CachedProperty.__get__
  File "/usr/local/lib/python2.7/site-packages/pandas/core/groupby/ops.py", line 252, in group_info
    comp_ids, obs_group_ids = self._get_compressed_labels()
  File "/usr/local/lib/python2.7/site-packages/pandas/core/groupby/ops.py", line 272, in _get_compressed_labels
    return compress_group_index(group_index, sort=self.sort)
  File "/usr/local/lib/python2.7/site-packages/pandas/core/sorting.py", line 379, in compress_group_index
    comp_ids, obs_group_ids = table.get_labels_groupby(group_index)
  File "pandas/_libs/hashtable_class_helper.pxi", line 1255, in pandas._libs.hashtable.Int64HashTable.get_labels_groupby
  File "pandas/_libs/hashtable_class_helper.pxi", line 199, in pandas._libs.hashtable.Int64Vector.resize
MemoryError: cannot allocate memory for array

The machine has 8Gb of RAM. Any ideas?

It seems to be some kind of a memory leak in pandas's read_csv, maybe due to it doing gzip-decompression.

Doing the decompression in a separate process and piping the result into cooler cload pairs should do the trick.

gzip -dc GSE109344_monkey_fibro_allValidPairs.txt.gz | cooler cload pairs -c1 2 -p1 3 -c2 5 -p2 6 /data/chromInfo.rheMac2.
simple.txt:${BINSIZE} - /data/${FILE}.cool


Running htseq-count over BAM files I am trying to derive an expression matrix from BAM files using htseq-count. These are bulk RNASeq BAM's by the way.
I have read the htseq-count documentation as well as samtools and figured that the following command should work:
samtools sort -on file.bam | samtools view -h - | htseq-count --type=exon --idattr=gene_id --additional-attr=gene_name - Homo_sapiens.GRCh38.77.gtf > counts.txt

The .gtf file here is from ensembl.
I have tried many variations of the above command but have failed to generate a count matrix. The error message states:
"0 SAM alignments processed"
Here are the first few lines of my .gtf file:
#!genome-build GRCh37.p13
#!genome-version GRCh37
#!genome-date 2009-02
#!genome-build-accession NCBI:GCA_000001405.14
#!genebuild-last-updated 2013-09
1       pseudogene      gene    11869   14412   .       +       .       gene_id "ENSG00000223972"; gene_name "DDX11L1";
gene_source "ensembl_havana"; gene_biotype "pseudogene";
1       processed_transcript    transcript      11869   14409   .       +       .       gene_id "ENSG00000223972"; trans
cript_id "ENST00000456328"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene"; transcript_nam
e "DDX11L1-002"; transcript_source "havana";
1       processed_transcript    exon    11869   12227   .       +       .       gene_id "ENSG00000223972"; transcript_id
 "ENST00000456328"; exon_number "1"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene"; trans
cript_name "DDX11L1-002"; transcript_source "havana"; exon_id "ENSE00002234944";

And all the genes in the text file have an expression value of 0. What exactly is going wrong in this step? Is the conversion to a SAM file the problem? I'm completely lost so any help would me much appreciated.
Thanks!
EDIT:: For anybody who sees this post cause they have a similar problem. Just note that the newer versions of htseq-count don't require sorted .bam files and, so following the editing of the .gtf file, all I needed to do was convert it to .sam using samtools view -h and then pipe this to htseq-count. Sorting the files prior to this conversion would cause errors.

If your BAM file is using UCSC chromosome names (chr1, chr2, etc.) and your GTF file is using Ensembl chromosome names (1, 2, etc.) then htseq-count won't know how to match your alignments to any of the genes. The simplest solution is to not use htseq-count, but instead to use featureCounts, which will handle the chromosome name difference transparently (it's one of the few programs that does) and is much faster anyway. In the future, download your GTF and fasta files from the same source, don't mix and match them.


qvalue and p.adjust functions Between the qvalue package and the p.adjust function, which is more appropriate to use when trying to calculate the q-values of a dataset? According to the manual for the q-value package, the q-value calculated is not an adjusted p-value which is what the p.adjust function would return. There seems to be an equal amount of people that use one or the other, but I wanted to see other opinions that could help clarify what would be most appropriate to use.

p.adjust can calculate adjusted p-values with different methods. I think you can choose Bonferroni correction and BH correction, for instance. BH correction gives you False Discovery Rate (FDR), which is one of the most commonly used I guess.

Q-value is different and cannot be calculated by p.adjust. The idea is that when you have many p-values to adjust, they form a distribution. If all p-values come from the null hypothesis, the p-value distribution is going to be uniform on [0, 1]. But usually, our p-values are a mixture of the null distribution and alternative distribution. The parameter $\pi_0$ describes the proportion of p-values (or tests) that comes from the null hypothesis. To calculated q-values, the algorithm first estimates $\pi_0$, then calculates the q-value.

Note that if you set $\pi_0=1$, which means that you believe all the p-values come from the null distribution, then q-value is equal to FDR.

I think most people do not know about the meaning and details of q-values and it might be sufficient to use BH procedure in most cases.


Detection of CNV(InDel) of intermediate size My impression is that small InDel (a couple of bp) is identified through cigar string in BAM and typical CNV (at least thousands of bp) is detected through read depth.

What about InDel or CNV with size these them, say hundreds of bp? These CNV is too long to be covered by read and too short to be detected by statistically comparing the difference of read depth.

Do we have consensus that this kind of CNV is hard to detect? Or is there a common strategy of detection for them? Do we generally think them having impact on protein function (especially in oncology background) ?

My impression is that small InDel (a couple of bp) is identified through cigar string in BAM...

Well, yes, by gaps in the alignment, easily detected by examining the CIGAR string.

  ...and typical CNV (at least thousands of bp) is detected through read depth.

Also by deviations from expected read pair distance or orientation.

  What about InDel or CNV with size these them, say hundreds of bp? These CNV is too long to be covered by read and too short to be detected by statistically comparing the difference of read depth. Do we have consensus that this kind of CNV is hard to detect?

Yes, what you describe has certainly been recognized and acknowledged in the literature. It was one of the primary motivations for a recently published paper I contributed to on de novo (germline) variant discovery. Here's an excerpt from that paper:

  In a reference-mapping context, calling indels with confidence requires accurate mapping of each read spanning the indel, with all gaps arranged consistently. This is possible only for short indels and tends to be prone to error and misalignment. Thus prediction of indels with length 10 bp has proved to be very challenging and accompanied by high false-positive and false-negative rates. Furthermore, the prediction of SVs via read mapping is only possible through indirect signatures such as alterations in read depth or read-pair signatures. These signatures can be quite noisy and result in high rate of false-negative and false-positive prediction.

As for your question:

  is there a common strategy of detection for them?

Numerous mapping-free k-mer based methods for variant detection are coming up in the literature, and all seem to perform pretty well. The trick is how to define which k-mers are interesting—then you grab those k-mers (or the reads containing them) and assemble them into variant-spanning contigs. A few examples:

In the Kevlar paper, we defined  "interesting" k-mers as those with high abundance in an individual of interest but zero (or very low) abundance in the individual's parents.
The HAWK paper, which describes a k-mer based GWAS-like method, defines "significant" k-mers as those significantly over- or under-represented in one class of samples versus a second class.
The NovoBreak paper compares paired tumor and normal WGS data to find k-mers that span mutatino breakpoints in cancer genomes.

  Do we generally think them having impact on protein function (especially in oncology background)?

There is definitely potential that indels can have a substantial impact on protein function. But since their detection has been so difficult for so long, a lot remains unknown about their actual collective impact on function.


How to proceed after contigs.fa file obtained from VELVET assembly I want to assemble a linear genome and I have a reference genome too.

I ran the following commands for de novo velvet assembly:

velveth auto 31,45,2 -fastq -shortPaired1 -separate 2-CA-SVA-2_S2_L001_R1_001.fastq 2-CA-SVA-2_S2_L001_R2_001.fastq

I chose auto_33 for further analysis

velvetg auto_33 -exp_cov auto

Now my output files inside auto_33 folder are:

contigs.fa
Graph2
LastGraph
Log
PreGraph
Roadmaps
Sequences
stats.txt

But I have no clue how to move ahead. Can I do anything else with velvet or do I have to switch to another program?

The contigs.fa file is a FASTA formatted file of the final contigs generated by velvet.  I would start assessing the quality of the assembly that was generated. Here are some questions to get you started.  Really the next steps depend on how you plan to use the assembly.

What is the distribution of contig lengths? Average and N50 contig length

How does the assembly compare to the reference?

How much of the reference does the assembly cover?

When you align the original reads back to the contigs how well are they covered?

Do any of the read alignments not make sense?
Are there sections of the contigs that are not covered?


What kind of analysis is practically done on GSE data files? I have a GSE data file in csv file format containing fields such as: ID, adj.P.Val, P.Value, t, B, logFC, Gene.symbol, Gene.title. In which adj.P.Val, P.Value, t, B, logFC fields being numeric. What are the factors I need to consider if I want to cluster the data only on logFC using K-Means clustering algorithm ? And first of all is it feasible to perform clustering on GSE data files ? If yes, what should be the approach ? If not, what different kinds of analysis can be performed on such kind of datasets ?

Following is snapshot of the data:

Have you ever used GEO2R? You can use GEO2R to analyze the data dividing it in subgroups (clusters) that interest you.


Concerning the current outbreak, using bioinformatics Dear Researches community

A question concerning the current outbreak, recently we attempting to contribute to the study of the Covid19, with a tight submission time ( 4 weeks).

Inquiring help because I'm not a virologist with an insufficient biological background. However, I assume able to do in silico analysis or ( dry lab).
Moreover, new sample sequences from local strain will be provided. Assembly and annotation could be done at the first stage. Then perhaps a comparative genome investigation between the old SARS and the Cuvid19.

Furthermore, The virus estimation length that causes all this is only around 29,829 base pairs of RNA. The genome contains at least 10 predicted open reading frames (ORFs): ORF1a,  ORF1b, S,  3,  4a, 4b,  5,  E,  M  and  N,  with sixteen predicted nonstructural proteins being encoded by ORF1a/b.
I am wondering if someone can assist by considering a particular section to focus on or a study suggests. Perhaps we can contribute at least 0.001 % to the communities.

I have been reading articles around, and seems are thousand of them every day.

so I have come up with starting list:

Genome assemblies, annotation
epitome studies
Here you predict the protein sequences based on the genome that will sequence, then there are a number of tools that analyze the immunogenicity of different proteins
That can help you identify parts of the COVID-19 proteins that are likely to be detected by the immune system and they, of course, can be targeted by vaccine developers.
comparison between old saras Vs newCoVID
The M protein is a good target. It is extremely abundant in the nucleocapsid and is the protein responsible for giving the virus it's the shape, Almost everyone is looking at the S protein, which is justified. But there is not as much focus on the M protein.

Your response is highly appreciated if you could attach the any reference related to the list above.

My personal approach to something like this is to first try to reproduce existing work in order to learn to use relevant tools and understand the concepts involved. Then as you get more comfortable with the subject matter, you may spot an opportunity to add value.

If you'd like to get started with the phylogenomics side of things, you can use this public workspace developed by the viral genomics group at the Broad Institute, which includes viral genome sequences from NCBI SRA as well as some workflows for doing genome assembly and quality control: https://app.terra.bio/#workspaces/pathogen-genomic-surveillance/COVID-19


Why does IGV turn a VCF QUAL of "." into -10? IGV is turning QUAL "." into -10 from this VCF of structural variations.  I'm wondering if this means anything or if I should report a bug to the IGV devs.

Here's a row from a VCF:

#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT SAMPLE

chrY    11292098    DUP00008569 C   <DUP>   .   PASS    IMPRECISE;SVTYPE=DUP;SVMETHOD=EMBL.DELLYv0.8.1;CHR2=chrY;END=56834684;PE=3;MAPQ=24;CT=5to3;CIPOS=-265,265;CIEND=-265,265    GT:GL:GQ:FT:RCL:RC:RCR:CN:DR:DV:RR:RV0/1:-5.71458,0,-20.0913:57:PASS:0:90774:0:-1:6:5:0:0

You see the QUAL value is "."

IGV reports it as -10

(This call is probably a false positive but it's just an example. There are others)

Looks like they are treating the "." as a value of 1.  According to the VCF spec the quality is "−10log10 prob".  Where are you getting your QUALs from?  I have never seen the quals as a "." and IGV probably needs to do some translation into an int but a QUAL of -10 makes no sense.

I would report it and let them explain.


Differential expression analysis when nested effects I have 3 tumour samples from 3 patients from an experiment, I also downloaded 10 normal samples from TCGA. My design is like this

> mycols
               condition batch
normal1           normal     1
normal2           normal     1
normal3           normal     1
normal4           normal     1
normal5           normal     1
normal6           normal     1
normal7           normal     1
normal8           normal     1
normal9           normal     1
normal10          normal     1
OESO_036_a_RNA    cancer     2
OESO_013_a_RNA    cancer     2
OESO_005_a_RNA    cancer     2
>

But DESeq2 returns error when trying for differential expression analysis

> dds <- DESeqDataSetFromMatrix(countData=counts, colData=mycols, design=~batch+condition)
Error in checkFullRank(modelMatrix) :
  the model matrix is not full rank, so the model cannot be fit as specified.
  One or more variables or interaction terms in the design formula are linear
  combinations of the others and must be removed.

  Please read the vignette section 'Model matrix not full rank':

  vignette('DESeq2')
In addition: Warning message:
In DESeqDataSet(se, design = design, ignoreRank) :
  some variables in design formula are characters, converting to factors

Then how I could extract deferentially expressed genes between tumour and normal samples while they are coming from two different batches?

Any help please?

llrs is right.  "Matrix Model not full rank" means that two of your conditions are confounded.  You can't get your normals from one place and your experimentals from another

And no, there generally is no clever way out of this.  For all you know, any differences you see between your normals and your cancers are 100% due to being sequenced somewhere else, and have nothing to do with real biology.

If you want to run this through DESeq, drop batch effect.  But tell whoever you are reporting the results to that batch is hopelessly confounded with cancer state, and you can't tell which differences are due to what.


PASA pipeline: compare experimental transcripts to the reference annotation I would like to ask if anyone has experience in running a subset of the PASA pipeline, in particular for the reconciliation of some experimental 'transcripts' with the reference annotation.

In more detail, I am working with RNA-seq data from D. melanogaster. I have reconstructed the 'transcripts' using Trinity. I have aligned these 'transcripts' to the reference genome using GMAP. Now I would like to match these 'experimental transcripts' with the reference annotation, to see how they compare.

I was wondering if I can run just a subset of the whole PASA pipeline, so basically skipping the step of alignment to the reference and using the alignment file that I generated externally.

Also, I would be glad to receive more general redirecting to any released software allowing the comparison of experimental transcripts to the reference annotation.

Cross-posted on the PASA Google Group pasapipeline-users

Cross-posted on biostars

Answer from the authors:

  PASA ends up using gff3 format instead of sam format for uploading the
  alignments.  If somebody has the gmap gff3 output, then can be done to
  just upload that directly using the custom alignment importer, but not
  with sam. They suggest just having PASA rerun GMAP as part of its
  regular routine.


Viral genome assembly using broad viral ngs pipeline? I am trying to assemble RNA virus genome using Broad Viral NGS pipeline BROAD VIRAL NGS PIPELINE. I am two questions :

1) As this pipeline requires unaligned bam format as input, how do I convert fastq to bam ?

2) Some of the tools like   assemble_trinity require Trimmomatic clipDb. How can I get this clipDb ?

Since you are already using the Broad Tools sets you can use Picard FastqToSam to make the conversion

As far a clipDb I am unfamiliar with that and a quick google search and look at the trimmomatic manual were unhelpful


determining doublets in single-cell RNA-seq Doublets are a known problem with scRNA-seq experiments, where 2 or more cells are sometimes captured instead. To determine their presence, there are studies that mix multiple species (such as human and mouse) or distinct cell types and then count how many cells contain genes from the very different transcriptomes. In those experiments, you can pick any cell and determine if it's a doublet based on a few key genes. In a typical single-cell experiment on one sample, you don't have the luxury of knowing what cells to expect and if a particular cell looks abnormal.

If you have a large fraction of doublets, you should be able to see a double peak if you plot the number of genes/transcripts per cell. However, if you have something like 1% doublets, you probably wouldn't notice anything. It seems that most people simply filter out top X% cells based on the number of genes/transcripts. How do you know what X should be? Can you actually estimate the fraction of doublets in a given experiment?

Update:

I found a nice figure from Stoeckius et al. to illustrate the issue:

Although the number of transcripts/UMIs seems to be roughly double for doublets/multiplets (median of ~1000 vs ~500), there is still high overlap.

Expected rates of doublets / duplets / multiplets

Fluidigm C1 doublet rate: around 1-5% depending on chip type used.
More information: Fluidigm white paper: Redesign of C1 Medium-Cell 96 IFCs Improves Single-Cell Capture Efficiency (PN 101-3328 B1)
10x Genomics Chromium 3'mRNA-Seq doublet rate: around 1% per 1000 cells captured

Rates in the 'Chromium Single Cell 3’ Reagent Kits v2 User Guide':

Multiplet Rate (%)  # of Cells Loaded   # of Cells Recovered
0.40%   870 500
0.80%   1700    1000
1.60%   3500    2000
2.30%   5300    3000
3.10%   7000    4000
3.90%   8700    5000
4.60%   10500   6000
5.40%   12200   7000
6.10%   14000   8000
6.90%   15700   9000
7.60%   17400   10000

Detection and identification of doublets and multiplets

Experimental

Visual inspection (e.g. Fluidigm C1)
FACS index sorting data can be used to identify doublets.
Mixed species experiment.

Illumina | BioRad White paper Single-Cell RNA Data Analysis Workflow
Zheng et al., 2017: Massively parallel digital transcriptional profiling of single cells

Cell "hashing" with oligo-tagged antibodies (Stoeckius et al., 2017): cells (barcodes) positive for two or more oligos are multiplets.

Computational

Mixed gene expression (e.g. lymphoid and myeloid, or T and B cell)
Ilicic et al. describe an approach for detecting low quality cells (incl. multiplets). They train a support vector machine (SVM) model on a dataset annotated by microscopy inspection. They find, among others, that

"wells containing multiple cells (multiples) show similar expression and noise patterns to broken cells", suggesting "that multiple cells contain a mixture of broken and high quality cells."
difference in proportion of duplicated reads
higher ERCC spike-in ratios

demuxlet (Kang et al., 2017, Nature Biotechnology)
DoubletDetection


How to retrieve SNPs data of different humans? How can I retrieve only SNPs from the dbSNP's FTP site. The files are gzipped and organized by human chromosome (i.e, chr[1-22,X,Y]), are these for one human? I need only SNPs data for various humans.

If you mean that you want snps from individuals, instead of all together, you can find them in 1000 genomes. Here different individuals from different populations are sequenced and variants are called, including snps.

Maybe you should add to your question what your goal is, what would you like to find?


Affymetrix tags (same ID's) present in different places of the genome could you explain me how to interpret such situation?

Same tag aligned in several distint places, near different genes...
Is the annotation able to distuinguish such 'same' tags, during DGE analysis?
Why tags arent unique?

Is the annotation able to distinguish such 'same' tags, during DGE analysis?

No, you'll end up discarding such probes (assuming you have a reason to actually use microarrays still).

  Why aren't tags unique?

Platforms such as HG-U133 date back to 2002 (meaning they were designed a few years before that), when the human reference genome (and the reference transcriptome) were not nearly as "finished" as they are today. Probes that may have made sense given what we knew almost 20 years ago may no longer makes sense now given what we know. Further, microarrays are rather conservative, in that they tend to primarily gain probes over time and not have them removed (lest companies lose customers who are relying on given probes).


What is the best way to get a large number of RNA seq data from SRA in Python without being denied access I have the following code to download data from the SRA using multithreading in Python. After running this a few times now (for testing purposes), I keep getting denied access to the data. Not sure how to fix this. In particular, cating the output files gives:
<?xml version="1.0" encoding="UTF-8"?><Error>AccessDenied<Message>Access Denied</Message><RequestId>cb0f0e98-cafb-1dd7-9b7b-d8c49756ec52</RequestId><HostId>QeDGVwBXYp61J0B4_OUTn7UsEsiQEec0n18DAeR0kaE</HostId></Error>

Here is my code:
import threading
import requests
import time

start = time.perf_counter()

class MyThread(threading.Thread):
    def __init__(self, url):
        threading.Thread.__init__(self)
        self.url = url
        self.result = None
        self.filename = url.split('/')[-1]
    def run(self):
        res = requests.get(self.url)
        with open(self.filename, 'wb') as f:
            f.write(res.content)
urls = [
        'https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos1/sra-pub-run-1/SRR000001/SRR000001.1',
        'https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos1/sra-pub-run-1/SRR000001/SRR000001.2',
        'https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos1/sra-pub-run-1/SRR000002/SRR000002.1',
        'https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos1/sra-pub-run-1/SRR000002/SRR000002.2']

threads = [MyThread(url, ) for url in urls]
for thread in threads:
    thread.start()
for thread in threads:
    thread.join()

finish = time.perf_counter()

print(f'Finished in {round(finish-start, 2)} second(s)')

I once needed to download more than 1000 files from SRA, and I found a blog post that is very useful (https://reneshbedre.github.io/blog/fqutil.html). The idea is to use aspera to fetch SRA files. This is fast becasue of aspera. Then use fasterq-dump to convert SRA file to fastq. fasterq-dump is faster than fastq-dump partly because it does not compress the output file. I think compressing is the rate-limiting step here. After this, you can use pigz, which is a multi-threading version of gzip to compress the fastq files. Note that if you want to download files to a computing cluster, you can submit jobs for fasterq-dump and pigz. This makes things much faster. You cannot submit jobs for fetching SRA, because you need to connect to internet.


Genomic coordinates of predicted binding pocket I would like to know the genomic coordinates of a proteins ligand-binding pocket.  For example, on this PDB page you can see that there is xray crystallography data for how CYP2C9 interacts with two ligands. I would like to know which amino acids form the space on the protein that binds the ligands, then translate (de-translate?) that back to the nucleotide sequence. Is anyone familiar with either a database that has this information or how to extract the ligand interaction data from the PDB files?

Thank you

It turns out there is a tool that does this beautifully.  I used the proteinToGenome function from the ensembldb package.

The code looks like this:

# Define the Ensemble database object for the
edbX <- filter(EnsDb.Hsapiens.v75, filter = ~ seq_name == "X")
# Set the first and last amino acid you want to convert and the ensembl protein id
GAGE10_prt <- IRanges(start = 5, end = 9, names = "ENSP00000385415")
# Perform the conversion.  GAGE10_prt is an object that has the coordinates
GAGE10_prt <- proteinToGenome(GAGE10_prt, edbX)


Unexpected error using "stdin" in bedtools intersect as a piped command From the bedtools intersect man (for version 2.30.0):
-a BAM/BED/GFF/VCF file “A”. Each feature in A is compared to B in search of overlaps.
   Use “stdin” if passing A with a UNIX pipe.

The stdin part, however, does not apply to version 2.27.1. The help on this version gives no advice for piping.
bam mem ref fq1.fastq fq2.fastq | samtools view -b | bedtools intersect -a stdin -b blacklist.bed -v > align_black.sam
Error: unable to open file or unable to determine types for file stdin

- Please ensure that your file is TAB delimited (e.g., cat -t FILE).
- Also ensure that your file has integer chromosome coordinates in the
  expected columns (e.g., cols 2 and 3 for BED).

I also tried:
(map to bam) | bedtools intersect -a "stdin" -b blacklist.bed -v > align_black.sam
(map to bam) | bedtools intersect -a -b blacklist.bed -v > align_black.sam
(map to bam) | bedtools intersect -b blacklist.bed -v > align_black.sam

But these don't work either, either for the same or a different reason. Here is a site that tells me to use the above syntax: https://bedtools.readthedocs.io/en/latest/content/example-usage.html
What is the solution to this if I want to pipe everything together for my users?

With v2.27.1 (or any version I ever used) instead of stdin use -
(...) | bedtools intersect -a - -b blacklist.bed -v > out.bam

Be aware that the output of this is a BAM, not SAM file. It would be easier though to simply make a list of regions you want to keep, e.g. the complement of the backlist file with the entire genome (bedtools complement), and then use samtools view -L option to only keep overlapping reads. That saves you the run through bedtools.


Java error when launching Pilon I am trying to use Pilon to improve a reference based on some Illumina data I have got. So, I aligned the Illumina reads to the reference using bwa. Then, I want to use Pilon to improve the reference using the alignment file in BAM format.

Information about Pilon is found here. I installed Pilon using conda.

The output of 'pilon' is:

Pilon version 1.22 Wed Mar 15 16:38:30 2017 -0400

    Usage: pilon --genome genome.fasta [--frags frags.bam] [--jumps jumps.bam] [--unpaired unpaired.bam]
                 [...other options...]
           pilon --help for option details

The output of:

pilon --genome ref.fa --frags readsSorted.bam --output pilo --outdir output/scaffold

is:

 Exception in thread "main" java.lang.reflect.InvocationTargetException
    at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
    at sun.reflect.NativeMethodAccessorImpl.invoke(NativeMethodAccessorImpl.java:62)
    at sun.reflect.DelegatingMethodAccessorImpl.invoke(DelegatingMethodAccessorImpl.java:43)
    at java.lang.reflect.Method.invoke(Method.java:498)
    at com.simontuffs.onejar.Boot.run(Boot.java:340)
    at com.simontuffs.onejar.Boot.main(Boot.java:166)
Caused by: java.lang.OutOfMemoryError: Java heap space
    at org.broadinstitute.pilon.GenomeRegion.<init>(GenomeRegion.scala:55)
    at org.broadinstitute.pilon.GenomeFile$$anonfun$contigRegions$1.apply(GenomeFile.scala:72)
 at org.broadinstitute.pilon.GenomeFile$$anonfun$contigRegions$1.apply(GenomeFile.scala:72)
    at scala.collection.TraversableLike$$anonfun$map$1.apply(TraversableLike.scala:234)
 at scala.collection.TraversableLike$$anonfun$map$1.apply(TraversableLike.scala:234)
    at scala.collection.immutable.Range.foreach(Range.scala:160)
    at scala.collection.TraversableLike$class.map(TraversableLike.scala:234)
 at scala.collection.AbstractTraversable.map(Traversable.scala:104)
 at org.broadinstitute.pilon.GenomeFile.contigRegions(GenomeFile.scala:72)
 at org.broadinstitute.pilon.GenomeFile$$anonfun$2.apply(GenomeFile.scala:52)
    at org.broadinstitute.pilon.GenomeFile$$anonfun$2.apply(GenomeFile.scala:52)
    at scala.collection.immutable.List.map(List.scala:277)
    at org.broadinstitute.pilon.GenomeFile.<init>(GenomeFile.scala:52)
    at org.broadinstitute.pilon.Pilon$.main(Pilon.scala:98)
    at org.broadinstitute.pilon.Pilon.main(Pilon.scala)
    ... 6 more

Like gringer already commented, the cause of your error is nicely given by the java program. Your default max heap settings are probably too low for the program. If you have more memory on your machine you can give max heap a higher value. See this post and answers on Stack overflow.

You'll probably need to use the -Xmxn argument in java.


processing fastq files using cellranger on linux I am trying to perform a cellranger count on fastq files generated from a 10x genomics single cell RNA Seq run.

Just to provide some background, I ssh’d into the AWS using our AWS IP with MobaXTerm linux command line.

My files are in the Isilon. The file folders are path/to/runs/ - this folder contains the .fastq files for each sample (ex: SI_GA_A1_S1, SI_GA_B1_S2, SI_GA_C1_S3, and SI_GA_D1_S4). I can see these folders listed in the linux command line.

The cellranger is found when I write the following:

path/to/cellranger-3.0.1/cellranger-3.0.1

then ls brings up 13 items located in the folder, one of them being cellranger.

When I write the following command,

/path/to/cellranger-3.0.1/cellranger-3.0.1/cellranger count \
    --id=run2\
    --transcriptome=/path/to/refseq/refdata-cellranger-hg19-3.0.0\
    --sample=sample1\
    --expect-cells=5000

and it gives me:

-bash: /path/to/cellranger-3.0.1/cellranger-3.0.1/cellranger: No such file or directory

This looks correct from what the 10x website says for single cell RNA Seq. Please advise.

$ ls -l /path/to/cellranger-3.0.1/cellranger-3.0.1/cellranger
-rwxrwx--- 1 1003080 1000000 460 Dec 5 19:41 /path/to/cellranger-3.0.1/cellranger-3.0.1/cellranger

this has been resolved. I downloaded the windows file and this caused trouble but when downloading the file within the Linux, it seems to be working.


How to merge .fastq.qz files into a single .fastq.gz with their same id without losing any content in parallel I have a large number of .fastq.gz files of different lane and reads. I have to merge them each reads group files into single .fastq.gz files.

**eg:

1st type
NA24694_GCCAAT_L001_R1_001.fastq.gz
NA24694_GCCAAT_L001_R1_002.fastq.gz
…....
…....
NA24694_GCCAAT_L001_R1_033.fastq.gz

2nd type
NA24694_GCCAAT_L001_R2_001.fastq.gz
NA24694_GCCAAT_L001_R2_002.fastq.gz
…....
…....
NA24694_GCCAAT_L001_R2_027.fastq.gz

3rd type
NA24694_GCCAAT_L002_R1_001.fastq.g
NA24694_GCCAAT_L002_R1_002.fastq.g
…....
…....
NA24694_GCCAAT_L002_R1_040.fastq.g

4th type
NA24694_GCCAAT_L002_R2_001.fastq.gz
NA24694_GCCAAT_L002_R2_002.fastq.gz
…....
…....
NA24694_GCCAAT_L002_R2_040.fastq.gz

so now i need to merge all the corresponding files into single file.

Output:
EA00694_GCCAAT_L001_R1.fastq.gz
EA00694_GCCAAT_L001_R2.fastq.gz
EA00694_GCCAAT_L002_R1.fastq.gz
EA00694_GCCAAT_L002_R2.fastq.gz

I tried ""cat'" but its not possible to do for large files. Could any one help me to run in script parallel in perl or python or shell. or any linux command which can be used for large file.

Thanks all.

All you need is cat. You won't find any better tool for a simple job like this. Just run:

cat NA24694_GCCAAT_L001_R1*fastq.gz > EA00694_GCCAAT_L001_R1.fastq.gz
cat NA24694_GCCAAT_L001_R2*fastq.gz > EA00694_GCCAAT_L001_R2.fastq.gz
cat NA24694_GCCAAT_L001_R3*fastq.gz > EA00694_GCCAAT_L002_R1.fastq.gz
cat NA24694_GCCAAT_L001_R4*fastq.gz > EA00694_GCCAAT_L002_R2.fastq.gz

That should work just fine. You say you tried cat and it "didn't work" but since you don't tell us how it failed, I can't really help. The only issues I can think of is that either there are too many files (as in several hundred thousand, whatever the value returned by getconf ARG_MAX on your system which is 2097152 on mine) or, more likely, you are running out of disk space.

If it's a disk space issue, you might be able to get around it by adding each file and then deleting it:

for i in 1 2 3 4; do
    for file in NA24694_GCCAAT_L001_R${i}_*fastq.gz; do
        cat "$file" >> EA00694_GCCAAT_L001_R${i}.fastq.gz && rm "$file"
    done
done

IMPORTANT: the command above is destructive. It will delete each file after it has been added to the new one. If everything works fine that's not a problem, but if I've made a mistake or your file names are slightly different you might lose your data. So I strongly urge you to simply run the cat commands in the first section on a machine with enough disk space to deal with them.


Can I eject the flow cell MinION from USB port when doing basecalling? Can I eject the flow cell MinION from USB port when doing basecalling? Because its taking long and its not possible to wait for its finished, will it affects basecalling?

You can, you will lose the transfer from the MinION. It would probably be better to just stop the run. The computer will continue to basecall untill it is done. All the raw reads are stored locally on the machine so it should be fine.


How to assess the quality of assembled .fasta genome files? I have assembled 3 .fasta files from contigs infastq format of 3 different Homo sapiens.

I would like to see if the assembled genomes which are in .fasta format are in good quality or not. In addition, I would be happy to see heritages of genes from the maternal or the parental genomes.

More specifically, I have 3 .fasta files of a mother, father and daughter.

How can I see they are indeed the biological parents?
How can I study what genes were inherited from whom of the parents?
How can I assess the quality of the assembled genomes? Is there a standard way to do that?
Is there a simple frequency analysis of bases n-grams I can use?

To expand on the comment by Devon Ryan, what you really need are (A.) bam alignment files for each individual and (B.) VCF files for each individual showing all of their genomic variants (with respect to the reference human genome).  You mentioned that you had fastq files and that you converted them into fasta files.  Do you know if the bam files and VCF files are also available for download anywhere?

If the bam files are not available you will have generate them using sequence aligner software (for example: bwa-mem).  With the bam files you can calculate the average depth of coverage (and many other things).

If the VCF files are not available then you will have to generate them from the bam files using variant calling software (for example, the GATK package from the Broad Institute).

In terms of assessing the quality of the sequence runs, a commonly applied tool is name FastQC, but once again, it uses fastq formatted files as input, not fasta files.

The main use for fasta formatted files these days seems to be for compact storage of various reference genomes from different genome sequencing projects (using gzip or similar type of compression).

With the VCF files for each individual in hand you should be able to begin addressing your other questions regarding the frequency of shared variants, etc.


Pathway enrichment analysis with or without custom background I am doing a pathway enrichment analysis of a sets of proteins.
Specific cell types (cardiospheres) were used to perform the experiment and for the quantification of proteins was used mass spectrometry.
I am using gProfiler to perform pathway analysis, should I use custom background of all of those proteins identified from the experiment? or should I used without custom background?
Thank you!
UPDATE: Also, I am doing an enrichment analysis of the interactome of the set of proteins of interest. In that case should I used a custom background for the enrichment analysis or not?
Thanks!!

I would defninitely go with a custom background which should be all the proteins that were actually measured in your experiment since you should only take into account what was measured.


Counting the characters in a row using python I am using python 3 and biopython 1.72. I have been using a protein sequence file in clustal format.

CLUSTAL 2.1 multiple sequence alignment

Human           ------------------------------------------------------------
Chimpanzee      ------------------------------------------------------------
Dog             ------------------------------------------------------------
Mouse           ------------------------------------------------------------
Xenopus         ------------------------------------------------------------
Amphioxus       MQWTGFRVSMTTLMMIMGVVAVLIALLPAKAQQPHDKSLRTTSTLTDTGASADEADMGSA

Human           ------------------------------------------------------------
Chimpanzee      ------------------------------------------------------------
Dog             ------------------------------------------------------------
Mouse           ------------------------------------------------------------
Xenopus         ------------------------------------------------------------
Amphioxus       HVELLDGDDDVGNGSDQMMVTLHLQSIFQCIRRPCEKVDRAIDPVTQRWRTANTRNDYQK

Human           ------------------MRLRVRLLKRTWPLEVPETEPTLGHLRSHLRQSLLCTWGYSS
Chimpanzee      ------------------MRLRVRLLKRTWPLEVPETEPTLGHLRSRLRQSLLCTWGYSS
Dog             ------------------MKLRVRLQKRTWPLDLPDAEPTLGQLRAHLSQALLPSWGFGS
Mouse           ------------------MKLRVRLQKRTQPLEVPESEPTLGQLRAHLSQVLLPTLGFSS
Xenopus         ------------------MKLRVRVRKQTNRLELEAESPTLGDLRSKLSSVTLPALGYST
Amphioxus       INVCVVPAYDVSLSTGVRMKLRVKISGQKTRVDVGQDCHTLGTLRTLLAPVLGEQYGLGD
                                  *:***::  :.  :::     *** **: *        * .

Human           NTRFTITLNYKDPLTGDEETLASYGIVSGDLICLILQDDIPA---------------PNI
Chimpanzee      NTRFTITLNYKDPLTGDEETLASYGIVSGDLICLILQDDIPA---------------PNI
Dog             DTRFAITLNNKDALTGDEETLASYGIVSGDLICLILEDAIPA---------------PNL
Mouse           DTRFAITLNNKDALTGDEETLASYGIVSGDLICLVLEDDMPA---------------PNL
Xenopus         EANFTITLNGKDALTGDQNTLESAGIISGDLIVVVLPDSQLP---------------PPQ
Amphioxus       DMPFEISLNGRDALLGDDKPLSDLGIVSGDLIHILLASVDQPSTNHNTQQQGQHPSSPEH
                :  * *:** :*.* **::.* . **:***** ::* .   .               *

I need to iterate over the rows of my file so that I could count the occurrence of characters in each row.
I tried using the code:

file_var = input("File: ")
count = 0
with open (file_var , 'r') as fo:
    for i in fo:
        count += 1
        print (i)
    print (count)

But this results just counting the total number of rows in a file(count variable). If I iterate over the file using biopython, it does'nt work also.

for seq_record in SeqIO.parse("sample.aln","clustal"):
    print (seq_record.id)
    print (repr(seq_record.seq))

seq_record.id results the names of the sequences used in the file.
repr(seq_record.seq) results the sequences of each specie alongwith description.

While I need to

get the rows only with amino acids (exclusing the symbols '*,.,:' line) and nothing else from the row and for each row I need to count the number of characters (exclusing '-' dashes)
the number of characters of that row be placed at the end of each of the row.

for example: for above example file, the counter should work like;

------------------------------------------------------------ 0
MQWTGFRVSMTTLMMIMGVVAVLIALLPAKAQQPHDKSLRTTSTLTDTGASADEADMGSA 60
------------------MRLRVRLLKRTWPLEVPETEPTLGHLRSHLRQSLLCTWGYSS 42


I think I might know what the problem is.

I plotted the log2 corrected values, and the distribution may not be normal because of the 0 values present in RNASeq data. If I understand correctly, this data is actually missing, whereas with my current analysis it is viewed by R as an actual value of 0. So either the missing data needs to be imputed, or the data needs to be scaled and normalized to mean of 0, then subsequently log2 transformed.

> data <- standardNormalization(data)
> data <- log2(data + 1)
> combat_data = ComBat(dat=data, batch=batch, mod=modcombat, par.prior=TRUE, prior.plots=FALSE)
> combat_data[1:5,1:2]
          X01240896.3f3f.4bf9.9799.55c87bfacf36
100507661                           -0.13857428
57103                               -0.06700506
22838                                0.16295096
55567                               -0.13579991
6147                                 1.64990772
          X01ad5016.f691.4bca.82a0.910429d8d25b
100507661                           -0.13789534
57103                               -0.06210232
22838                                0.01546224
55567                               -0.13443289
6147                                 1.48473072


Why don't the "child" peptide masses sum to the parent in this De novo sequencing example? In that example, 129+216+327 does not equal the parent's mass of 492.
My understanding is that during MS, the protein is broken down into peptide masses of varying length of approx. 5 to 20 amino acids. They also lose ions such as water or ammonia like in the example., the experimental spectra I thought would show the masses of these "children" that can be summed up to the parent (if you also add all the ammonia and water ions that are broken off). But that is not the case. Where am I wrong in my understanding?

In a peptide tandem mass spectrometry experiment, each fragment (child) ion is typically the result of breaking a single bond in the peptide backbone yield b and y ions. Because each ion fragment ion is the result of one fragmentation event, we don't expect the fragment ion masses to sum to the precursor mass.
If your peptide was fragmented well, you should get a ladder of b and y ions where the differences in masses between ions in the b and y ion series correspond to the masses of amino acids (+/- any modifications and neutral losses, such as the NH3 and H2O mentioned above).
A great resource for understanding these concepts is this tutorial:
https://www.ionsource.com/tutorial/DeNovo/introduction.htm


How to output all sequences with bwa mem, not `*`? I've been running bwa mem -a for alignment, using the -a flag---this will

  output all alignments for SE or unpaired PE

I've noticed in the SAM that there are several alignments with * in the SEQ and QUAL fields. Based on the documentation:

  SEQ: segment SEQuence. This field can be a ‘*’ when the sequence is not stored. If not a ‘*’, the length of the sequence must equal the
  sum of lengths of M/I/S/=/X operations in CIGAR. An ‘=’ denotes the
  base is identical to the reference base. No assumptions can be made on
  the letter cases.
  QUAL: ASCII of base QUALity plus 33 (same as the quality string in the Sanger FASTQ format). A base quality is the phred-scaled base
  error probability which equals −10 log10 Pr{base is wrong}. This field
  can be a ‘*’ when quality is not stored. If not a ‘*’, SEQ must not be
  a ‘*’ and the length of the quality string ought to equal the length
  of SEQ.

it appears the sequence isn't stored.

I would strongly prefer to have the sequence in this case. Is there any to direct bwa to output these sequences?

I don't believe that is an option.  Internally we created a patch and submitted it to BWA however that was 4 years ago and it has not been accepted.


Analyzing Illumina Counts I'm pretty new to all of this--forgive me if this is a simple question.

When I download illumina counts from GEO (like the supplementary file in
GSE89225). Can I do comparisons directly on that file? Is there some normalization procedure I should go through?

Thanks!

The counts files for GSE89225 is the output of HTSeq-count as a large matrix. Unless you are developing a differential expression package yourself you should not attempt to directly use this. Rather, you should load it into R and use packages such as DESeq2, edgeR, or limma (those are the most popular ones).

For convenience, in DESeq2 you would want the DESeqDatasetFromMatrix() function after loading this.


How can I edit a specific FASTQ read in place, given the read ID? Given a read ID, I want to edit a single basecall (e.g. the 12th base) for just that read from within a large FASTQ file containing millions of reads.

example, i want to change the 12th base ('C') in read 31027 to a 'T':

@70630 1:N:0:ATCACG
GAAGGTCCATGGATAATACTCAATTTTCCACAACAGCTTTTGTACTCTAGATCATTGATATTTACCAAAAGTCACTCAAACTCATCCTATGCATAATTCTAGTCCACCAATCATGATATGATGGAGAACATGGTTGTAATCAGGAAGACAG
+
DDDDDHIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIHIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIGHHIIIIIIIIIIIIIIIIIIIIIIHIIHIHIIIHHIGIIIIHHHHIIH@
@31027 1:N:0:ATCACG
CAAAAGTCACTCAAACTCATCCTATGCATAATTCTAGTCCACCAATCATGATATGATGGAGAACATGGTTGTAATCAGGAAGACAGATAAAGCAGCAGACCAAAAGTAATCTGAGAAATTATATTTGAATCACTCAGATATACATCAAATA
+
DDBDDHC@GHFHHIIE@CEHHIHEFFCGFCH?HHEGFCCEHHCGH@HFHHHCHHIIIHEHIIHII@CECHIHIIIEECCDGEHFFHHHHEHHFHHGGIHHHDGHFHIIIGHHHHHHEHHIIHIIIGGHHHCHHHCHIIIHHIEH@GHIHIC
@87319 1:N:0:ATCACG
CAATTAAGCTTTGGCAACGGTGGTCAAGATGAGATGCATATGGAGATAATAACTAAAAGTCAATCGAGACTCATCGTATGCATATTTCTAGTCCATCGATCATGAAATGATAGGATAGCTAGAATGAAAAGTAAATTTCCAGAAGGTCCAT
+
DDDDDIIIIIIIIIIIIIHIHHIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIHHIIIIIIIIIIIIIIIIIHIIHIIIIIIIIIIIIHIIIIIIIIIHIHIH

Of course I could stream through the entire FASTQ file rewriting everything to a new file until i get to the read of interest, do the edits to that read, write it out to the new file, and then keep going until the end. BioPython would work fine for this as would other approaches.

However, is there an efficient way to do this without the full read/write stream? Could this be done with sed/awk ? What about an indexed fastq file?

cheers
David

You'll first need to determine the appropriate line number, which you can do with grep -m1 -nw "@31027" foo.fastq. After that, note that you can provide a line number to sed:

sed -i '123456s/CAAAAGTCACTCA/CAAAAGTCACTTA/' foo.fastq

That will do the replacement only on line 123456 and edit the file in place (the -i option).


Replace a prefix of a string I have a data frame like this

> head(a[1:2,1:4])
   Tumor_Sample_Barcode Chromosome Position End_Position
1:                   A1       chr4 90169866     90169866
2:                   A1      chr11 60235747     60235747
    > dim(a)
[1] 753655    134
>
> unique(a$Tumor_Sample_Barcode)
 [1] "A1"  "A10" "A11" "A12" "A13" "A14" "A15" "A16" "A17" "A18" "A19" "A2"  "A20" "A21" "A22" "A23" "A24" "A25" "A26" "A27" "A28"
[22] "A29" "A3"  "A30" "A31" "A32" "A33" "A4"  "A5"  "A6"  "A7"  "A8"  "A9"

I want to change a$Tumor_Sample_Barcode  from A1:A21 to res1:res21 but I don't how to do that

I think you are looking for gsub function:

a$Tumor_Sample_Barcode <- gsub("A","res",a$Tumor_Sample_Barcode)

If you prefer, you can use the str_replace_all function from stringr package. But it tooks a little bit more time than gsub:

Dummy data

seq = as.factor(paste0("A",1:21))
library(data.table)
dt <- data.table(Barcode = rep(seq,10000),
                 Value = rnorm(210000))

Test of gsub vs str_replace_all

> system.time(dt$Barcode <- gsub("A","res",dt$Barcode))
   user  system elapsed
  0.077   0.000   0.077

> system.time(dt$Barcode <- str_replace_all(dt$Barcode,"A","res"))
   user  system elapsed
  0.123   0.004   0.126
```


BLASTn query coverage discrepancy Using the BLASTn webtool with,
Query Seq: CAATTGCGAAGCCACATTATT
Subject Seq: CP031332.1
BLASTn returns a single hit to CP031332.1 with query coverage at 100%.
However if one looks at the alignment, the query coverage is clearly 17/21 ≈ 81%. Could someone explain to me the reason for this discrepancy?
Edit,
Web version: https://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastn&PAGE_TYPE=BlastSearch&LINK_LOC=blasthome
Output: https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Get&RID=RCK3F55Y114

To expand on the excellent answer by haci:
BLAST reports (a) the summary of alignments and (b) individual alignments, see below. Note that your link to the results expired (as NCBI BLAST results links tend to do). It is considered best practice to post the text version of the relevant portion of the results in the question itself.
The summary provides an overview of all the query sequence individual hits in the subject sequence. Because the individual hits (alignments, or ranges) together cover 100% of the query and each of these has the max identity = 100%, that is what it reports in the summary. It is not very intuitive for many users.
You should always look at the individual hits in detail. The above 100% coverage in the summary consists of ranges 1 and 14, shown below. The top hit (range 1) covers only 17 out of 21 nt with 100% identity, and this is the most relevant information from BLAST. The other hit (range 14) covers 10 out of 21 nt with 100% identity. This is the cause of the discrepancy.
These 2 alignments are not close together in the subject sequence. Depending on your goal, choose one or more of the individual alignments and consider alignment length and % identity of those hits individually, rather than summarized together as in the BLAST report summary. For most applications, the top hit (range 1) would be the most relevant alignment.
Program: BLASTN
Subject:Neisseria meningitidis strain M22814 chromosome, complete genome ID: CP031332.1(nucleic acid) Length: 2190201
Query #1: seq1 Query ID: lcl|Query_56011 Length: 21

Sequences producing significant alignments:
                                                                  Max    Total Query   E   Per.
Description                                                       Score  Score cover Value Ident  Accession
Neisseria meningitidis strain M22814 chromosome, complete genome  34.2   57729 100%  0.001 100.00 CP031332.1

Alignments:
>Neisseria meningitidis strain M22814 chromosome, complete genome
Sequence ID: CP031332.1 Length: 2190201
Range 1: 2074883 to 2074899

Score:34.2 bits(17), Expect:0.001,
Identities:17/17(100%),  Gaps:0/17(0%), Strand: Plus/Plus

Query  1        CAATTGCGAAGCCACAT  17
                |||||||||||||||||
Sbjct  2074883  CAATTGCGAAGCCACAT  2074899

...

Range 14: 354292 to 354301

Score:20.3 bits(10), Expect:15,
Identities:10/10(100%),  Gaps:0/10(0%), Strand: Plus/Minus

Query  12      CCACATTATT  21
               ||||||||||
Sbjct  354301  CCACATTATT  354292


Convert a fasta file to Newick format My problem is (should be) pretty simple.
I have a series of 10 sequences which I wish to convert into a Newick format file - what would be good software to do this? I've tried using Jalview and am getting nowhere.
Any help appreciated.

The best would be Biopython using AlignIO, given they are all in a single file otherise SeqIO. The alignment approach,
from Bio import AlignIO
align = AlignIO.read("/home/path/myseq.fa", "fasta")
print(align)
with open("/home/path/example.nex", "w") as handle:
    SeqIO.write(align, handle, "nexus")

(not debugged)
MEGAX or ClustalX which are nice GUI's that will also work. If you've not got Biopython installed then ClustalX is probably the easiest.


How to use C++ htslib to read VCF contig name and size? A typical VCF file has:

##contig=<ID=chr1,length=248956422>
##contig=<ID=chr2,length=242193529>

I would like to use htslib in C++ to read it. My attempt:

htsFile *fp = bcf_open("my.vcf", "r");
bcf_hdr_t *hdr = bcf_hdr_read(fp);

In https://github.com/samtools/htslib/blob/develop/htslib/vcf.h, I'm not able to find a function that can do that for me.

How to read chr1 and 248956422 in C++?

#include "htslib/vcf.h"

int main(int argc, char *argv[])
{
    htsFile *fp;
    bcf_hdr_t *hdr;
    bcf_idpair_t *ctg;
    int i;
    if (argc == 1) {
        fprintf(stderr, "Usage: print-ctg <in.vcf>\n");
        return 1;
    }
    fp = vcf_open(argv[1], "r");
    hdr = vcf_hdr_read(fp);
    ctg = hdr->id[BCF_DT_CTG];
    for (i = 0; i < hdr->n[BCF_DT_CTG]; ++i)
        printf("%s\t%d\n", ctg[i].key, ctg[i].val->info[0]);
    bcf_hdr_destroy(hdr);
    vcf_close(fp);
    return 0;
}

On stability: this use has been in htslib forever. In general, functions/structs/variables in the public headers are meant to be stable. However, there is no guarantee that future versions will always keep the same APIs.


How to convert a binary matrix of gene presence or absence into a fasta sequence I have a binary matrix of gene presence or absence which looks like: [roary output]

Gene    sample1 sample2 sample3 sample4
fliI        1       1       1        1
patB_1      1       1       1        1
pgpA        1       1       1        1
osmB        1       1       1        1
cspA        1       0       1        1

How can I convert to fasta so it looks like and technically aligned:

>sample1
11111
>sample2
11110

Is it possible to do this in R?

It's not a fasta file, but:

> m
       sample1 sample2 sample3 sample4
fliI         1       1       1       1
patB_1       1       1       1       1
pgpA         1       1       1       1
osmB         1       1       1       1
cspA         1       0       1       1
> # Collapse to labeled strings
> blah = apply(m, 2, function(x) paste(x, collapse=''))
> blah
sample1 sample2 sample3 sample4
"11111" "11110" "11111" "11111"
> # Write that to a file with the appropriate format
> cat(paste(mapply(function(x, y) sprintf(">%s\n%s\n", x, y), names(blah), blah), collapse=""), file="some file.fa")

Make sure to change some file.fa. I have no idea why you would want to do all of this, but this will give you the output you want:

$ cat blah.txt
>sample1
11111
>sample2
11110
>sample3
11111
>sample4
11111


How to find the number of contigs produced? N50 length in base pairs? N90 length in base pairs? I am having trouble with some underlying questions about my project.
I have ran a Trinity tool to create a trinity.fasta file.
To determine some underlying questions, I used the utility asm_stats to gather some basic statistics.
I am just having some issues of calculating data.
DATA:
Main genome scaffold N/L50: 7464/1.032 KB
Main genome contig N/L50: 7464/1.032 KB
Main genome scaffold N/L90: 14314/106
Main genome contig N/L90: 14314/106
Number of contigs (scaffold length):
25,534 (ALL)
25,534 (100)
21,609 (250)
10,818 (500)
907 (1 KB)
38 (2.5KB)
5 (5 KB)
2 (10KB)
With the data attached above, would my N50 and N90 length in BASE PAIRS so excluding the kb, be 7464 and 14314, respectively?
Would the number of contigs produced by this assembly be 7, since there are 7 different contigs produced or would it be 25,534?
Thanks!

With the data attached above, would my N50 and N90 length in BASE PAIRS so excluding the kb, be 7464 and 14314, respectively?

Looks like this is inverse so your N50 is 1032bp and your N90 is 106bp. N90 cannot be longer than N50 by definition. You L50 and L90 are 7464 and 14314

Would the number of contigs produced by this assembly be 7, since there are 7 different contigs produced or would it be 25,534?

25534
Now for the real answer. Do not calculate or report N50 on transcript assemblies. It is more than useless. Transcriptomes are - by nature - dynamic in length. Longer is not better in this regard. Are your short transcripts lower quality assemblies, probably not, but N50 treats them as such. Please read the Trinity wiki for other, better ways to QC your transcript assembly.


How to map PDB chains to Uniprot IDs using API services I have a lot of PDB IDs and I need to get uniprot fasta sequences of these PDB IDs special chains by API services. For example, imagine that I need to get fasta sequence of '1kf6' 'A' chain. The uniprot entry (Accession) for this '1kf6' chain is 'P00363'. If I find a way to get this entry, I can simply retrieve its fasta sequence using bioservices package:

from bioservices import *
u = UniProt()
sequence = u.retrieve("P00363","fasta")
print(sequence)

The main problem is that I don't know how can I find this uniprot entry("P00363" here) for each chains in PDB IDs(here I imagine that I don't know this "P00363" is the uniprot entry for '1kf6' 'A' chain and I'm trying to find it using API services). I tried to use ID mapping. In ID mapping cross-references, the chain name for each entry is specified and if I can add a cross-reference column to u.search('1kf6', columns="id,entry name,length, genes,pdb"), I can simply find the entry related to my favourite chain. The following code returns some entries each related to one of '1kf6' chains. The problem is that I don't know exactly which entry is related to 'A' chain that here is my favorite chain. In the following code, id searches for entry and other parameters search for sth else related to them:

from bioservices import *
u = UniProt()
res = u.search("1kf6",columns="id,entry name,length, genes")
print(res)

Each one of the following words can be used and added to columns part:

'id', 'entry name', 'genes', 'genes(PREFERRED)', 'genes(ALTERNATIVE)', 'genes(OLN)', 'genes(ORF)', 'organism', 'organism-id', 'protein names', 'proteome', 'lineage(ALL)', 'lineage-id', 'virus hosts', 'fragement', 'sequence', 'length', 'mass', 'encodedon', 'comment(ALTERNATIVE PRODUCTS)', 'comment(ERRONEOUS GENE MODEL PREDICTION)', 'comment(ERRONEOUS INITIATION)', 'comment(ERRONEOUS TERMINATION)', 'comment(ERRONEOUS TRANSLATION)', 'comment(FRAMESHIFT)', 'comment(MASS SPECTROMETRY)', 'comment(POLYMORPHISM)', 'comment(RNA EDITING)', 'comment(SEQUENCE CAUTION)', 'feature(ALTERNATIVE SEQUENCE)', 'feature(NATURAL VARIANT)', 'feature(NON ADJACENT RESIDUES)', 'feature(NON STANDARD RESIDUE)', 'feature(NON TERMINAL RESIDUE)', 'feature(SEQUENCE CONFLICT)', 'feature(SEQUENCE UNCERTAINTY)', 'version(sequence)', 'domains', 'domain', 'comment(DOMAIN)', 'comment(SIMILARITY)', 'feature(COILED COIL)', 'feature(COMPOSITIONAL BIAS)', 'feature(DOMAIN EXTENT)', 'feature(MOTIF)', 'feature(REGION)', 'feature(REPEAT)', 'feature(ZINC FINGER)', 'ec', 'comment(ABSORPTION)', 'comment(CATALYTIC ACTIVITY)', 'comment(COFACTOR)', 'comment(ENZYME REGULATION)', 'comment(FUNCTION)', 'comment(KINETICS)', 'comment(PATHWAY)', 'comment(REDOX POTENTIAL)', 'comment(TEMPERATURE DEPENDENCE)', 'comment(PH DEPENDENCE)', 'feature(ACTIVE SITE)', 'feature(BINDING SITE)', 'feature(DNA BINDING)', 'feature(METAL BINDING)', 'feature(NP BIND)', 'feature(SITE)', 'go', 'go(biological process)', 'go(molecular function)', 'go(cellular component)', 'go-id', 'interpro', 'interactor', 'comment(SUBUNIT)', 'citation', 'citationmapping', 'created', 'last-modified', 'sequence-modified', 'version(entry)', '3d', 'feature(BETA STRAND)', 'feature(HELIX)', 'feature(TURN)', 'comment(SUBCELLULAR LOCATION)', 'feature(INTRAMEMBRANE)', 'feature(TOPOLOGICAL DOMAIN)', 'feature(TRANSMEMBRANE)', 'annotation score', 'score', 'features', 'comment(CAUTION)', 'comment(TISSUE SPECIFICITY)', 'comment(GENERAL)', 'keywords', 'context', 'existence', 'tools', 'reviewed', 'feature', 'families', 'subcellular locations', 'taxonomy', 'version', 'clusters', 'comments', 'database', 'keyword-id', 'pathway', 'score'

I used some of the above words in columns part but they couldn't retrieve chain name. Anyone know which one of the above parameters can give me the chain name for this PDB ID? or anyone has any other idea for finding uniprot entry(UniprotKB AC) related to each PDB ID chain by API services?

I found the answer and I thought it may be useful for others. we can use pypdb package for that purpose (it works in unix systems):

import pypdb

all_info = pypdb.get_all_info('1kf6')
print(all_info)

the output is as follows:

{'polymer': [{'@entityNr': '1', '@length': '602', '@type': 'protein', '@weight': '66057.6', 'chain': [{'@id': 'A'}, {'@id': 'M'}], 'Taxonomy': {'@name': 'Escherichia coli', '@id': '562'}, 'macroMolecule': {'@name': 'Fumarate reductase flavoprotein subunit', 'accession': {'@id': 'P00363'}}, 'polymerDescription': {'@description': 'FUMARATE REDUCTASE FLAVOPROTEIN'}, 'enzClass': {'@ec': '1.3.5.4'}}, {'@entityNr': '2', '@length': '243', '@type': 'protein', '@weight': '27021.9', 'chain': [{'@id': 'B'}, {'@id': 'N'}], 'Taxonomy': {'@name': 'Escherichia coli', '@id': '562'}, 'macroMolecule': {'@name': 'Fumarate reductase iron-sulfur subunit', 'accession': {'@id': 'P0AC47'}}, 'polymerDescription': {'@description': 'FUMARATE REDUCTASE IRON-SULFUR PROTEIN'}, 'enzClass': {'@ec': '1.3.5.1'}}, {'@entityNr': '3', '@length': '130', '@type': 'protein', '@weight': '14898.8', 'chain': [{'@id': 'C'}, {'@id': 'O'}], 'Taxonomy': {'@name': 'Escherichia coli', '@id': '562'}, 'macroMolecule': {'@name': 'Fumarate reductase subunit C', 'accession': {'@id': 'P0A8Q0'}}, 'polymerDescription': {'@description': 'FUMARATE REDUCTASE 15 KDA HYDROPHOBIC PROTEIN'}}, {'@entityNr': '4', '@length': '119', '@type': 'protein', '@weight': '13118.9', 'chain': [{'@id': 'D'}, {'@id': 'P'}], 'Taxonomy': {'@name': 'Escherichia coli', '@id': '562'}, 'macroMolecule': {'@name': 'Fumarate reductase subunit D', 'accession': {'@id': 'P0A8Q3'}}, 'polymerDescription': {'@description': 'FUMARATE REDUCTASE 13 KDA HYDROPHOBIC PROTEIN'}}], 'id': '1KF6'}

This output gives all the information related to all chains and we can simply find uniprotkb AC of each chain in this output. I also found another answer for this question. This answer can be found here and the following code:

import requests
from xml.etree.ElementTree import fromstring

pdb_id = '1kf6.A'
pdb_mapping_url = 'http://www.rcsb.org/pdb/rest/das/pdb_uniprot_mapping/alignment'
uniprot_url = 'http://www.uniprot.org/uniprot/{}.xml'

def get_uniprot_accession_id(response_xml):
    root = fromstring(response_xml)
    return next(
        el for el in root.getchildren()[0].getchildren()
        if el.attrib['dbSource'] == 'UniProt'
    ).attrib['dbAccessionId']

def get_uniprot_protein_name(uniport_id):
    uinprot_response = requests.get(
        uniprot_url.format(uniport_id)
    ).text
    return fromstring(uinprot_response).find('.//{http://uniprot.org/uniprot}recommendedName/{http://uniprot.org/uniprot}fullName'
    ).text

def map_pdb_to_uniprot(pdb_id):
    pdb_mapping_response = requests.get(
        pdb_mapping_url, params={'query': pdb_id}
    ).text
    uniprot_id = get_uniprot_accession_id(pdb_mapping_response)
    uniprot_name = get_uniprot_protein_name(uniprot_id)
    return {
        'pdb_id': pdb_id,
        'uniprot_id': uniprot_id,
        'uniprot_name': uniprot_name
    }

print (map_pdb_to_uniprot(pdb_id))

The output of this code is as follows:

{'pdb_id': '1kf6.A', 'uniprot_id': 'P00363', 'uniprot_name': 'Fumarate reductase flavoprotein subunit'}


Confirmation about PERMANOVA and Variance partitioning (permutest) Could you please help me to understand what is the difference between Permutational multivariate analysis of variance (PERMANOVA,adonis) and  Variance partitioning
(permutest).
I understood that permutest is used to confirm the PERMANOVA result that differences in PERMANOVA is not due to dispersion.
###adonis_calculation
BC <- phyloseq::distance(phyloseq_prop, "bray")
pmv =adonis(BC ~Season*Region, as(sample_data(physeqN3), "data.frame"),permutations = 5000)
pmv

Call:
adonis(formula = BC ~ Season * Region, data = as(sample_data(physeqN3),      "data.frame"), permutations = 5000)

Permutation: free
Number of permutations: 5000

Terms added sequenti
Total         99   19.8578                 1.00000 ally (first to last)

              Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
Season         1    1.8586 1.85862  12.036 0.09360  2e-04 ***
Region         1    2.4218 2.42183  15.682 0.12196  2e-04 ***
Season:Region  1    0.7522 0.75222   4.871 0.03788  2e-04 ***
Residuals     96   14.8251 0.15443         0.74656

Many thanks

Permutational Multivariate Analysis of Variance (PERMANOVA) assumes no distribution which means it is a non-parameteric statistical test. According to the assumption of PERMANOVA test, each group must have same dispersion in the data. Before performing any statistical test, we must check the underlying assumptions. Variance partitioning (permutest) is a confirmatory test to check the assumption (similar dispersion). You can use permutest.betadisper for permutest. After that you can perform PERMANOVA test using adonis().
Interpretation of PERMANOVA test –

If p value is small (<0.05 ) we reject the null hypothesis (No difference between the groups) otherwise we accept the null hypothesis (Significant difference between the groups).

Coefficient of determination, R-squared value indicates the variation in the distances being explained by groups under the study. Higher value indicates a better fit.

You can further look into the confidence interval for each group.
Hope this confirms the answer to your question.


How can I output an identity matrix in progressiveMauve? I'm just getting started with the Mauve aligner and I'm finding the documentation a bit lacking. I'm using the progressiveMauve tool from the command line and would like to output an identity matrix file. The set of output files it generates by default for progressiveMauve --output=outfile are outfile, outfile.backbone and output.bbcols.

The documentation for the tool lists an --input-id-matrix option but this doesn't seem to have any effect. The documentation page for output files does describe an identity matrix file but doesn't say anything about how to generate it.

mauve documentation have a similar command but the usage is --id-matrix=<file>. So very similar to the progressivemauve command (just doesn't include --input). Try it and see if it works.

I think it's a bug for sure as several other people on biostar have complained.


Quantify low complexity regions in DNA sequences I have fasta files with multiple sequences. They are reads that mapped to the genome outside of  probe-targeted regions.
From a quick perusal, they appear to be repetitive and have low complexity. Is there a way I can quantify their complexity and use that as an explanation as to why those regions were mapped to?

"Complexity" is also known as entropy. A quick and dirty way is to use a file compression utility - zip will get good compression for a repetitive sequence.
You can also look at the k-mer distribution for some suitable k: in a repetitive sequence it will be biased. Here is a tool that give a local plot: https://rdrr.io/github/vsbuffalo/qrqc/man/kmerEntropyPlot.html


How to count the kmer occurrence in FASTA file considering overlapping and reverse complement? I am using count_overlap() for counting the kmers from Biopython .
Does it take the reverse complement of kmer into account?
I need to count reverse complement too

There's a lot detail missing here, but at a guess could you try,
from Bio import SeqIO
from Bio.Seq import Seq
for record in SeqIO.parse("ls_orchid.gbk", "genbank"):
    reversec = record.seq.reverse_complement())
    revcount = reversec.count_overlap('AT')
    print (revcount + count_overlap('AT'))

Essentially count_overlap only counts the forward strand if this is the Biopython method. Not debugged the above.


ATAC-seq peak annotation and downstream analysis I have 7 different cell types, 4 normal cell types and 3 disease cell types.
I have done until peak calling. Now I want to do peak annotation.

The sample I have is all biological replicates. Should I merge the peak file and then do peak annotation, or should I be doing individual peak annotation?

The way we do this (and I think its quite common?), is to merge all the files of the same condition (taking the same number of reads from each), and then calling peaks on the merged sample - so you would have two peak sets - one from normal and one from disease. You then do the union of the two peak sets.

This makes annotation much easier.

For downstream quantitative analysis, you can then quantitate each individual sample against the merged peak set.


Seurat for clustering bulk RNA-seq? Is it ever ok to use Seurat for clustering bulk samples?

I am looking at FPKM data from ~750 bulk RNA-seq samples generated using Cufflinks. As suggested for FPKM data, I manually input log transformed data to the @data slot [cd138_bm@data <- log(cd138_bm@raw.data + 1)] and skip the NormalizeData() function. I then use functions FindVariableGenes, ScaleData, RunPCA, FindClusters, RunTSNE, FindAllMarkers in their usual ways to find clusters & cluster markers. My clustering results are quite reasonable, and reflect published work clustering similar samples.

What are the potential pitfalls of using these Seurat functions on bulk data? In FindAllMarkers, would you recommend I use the "negbinom" test? (currently using wilcox) Any other arguments you would recommend changing from the default?

I'm not sure Seurat is the best tool for this as it was developed for single cell RNA seq data and there are a few intricacies of that type of data that are very different from bulk RNA seq.
For bulk there are really good packages available and corresponding workflows, e.g. limma, edgeR and DESeq2.

The main problems with Seurat for bulk RNA-seq:

Seurat expects counts as input - FPKM are not counts nor are they log counts or log norm counts. It's a specific normalisation method that takes into account gene length and library size and breaks the link between gene counts and variance.
Because of the distribution of reads from a single cell that you get from the sequencer Seurat scales all unique reads (i.e. UMIs) to 10,000 per cell and usually regresses out the nUMI count. Again, becuase you've got FPKMs this makes no sense.
The other functions you can use on bulk data, i.e. tSNE, finding clusters and markers (which is essentially DGE analysis). The one thing to note for the latter is that again because you've got FPKMs they don't follow the negative binomial distribution. Your best bet is limma trend but I'm not sure it's included in Seurat.


Help interpreting my phylogenetic tree construction of bacterial species in same genus I constructed a phylogenetic tree using mole blast, and it created a neighbor joining tree, so I used the alignment of those sequences with mine (1-43) in mega with TN93+G+I model but I don't understand how to interpret why my species in question (29) is so much different (or maybe novel?).

Attached is my original tree (Picture 1) and my bootstrap analysis (Picture 2):

The top phylogeny is likely mid-point rooted, but I suspect the rooting is fine. The lower diagram has been rooted on the query sequence. It is taxonomic analysis using 16S which in bacteria is considered to represent species delineation at a given percentage of divergence. I thought it was 5%, which is impossible based on the phylogram in the top picture.

Bootstrap values lower than 80% mean that the clade is not supported (monophyletic group is a better term). The delineation between 10001 and 10013 against the massive ingroup of 13 16S sequences is robustly supported by bootstraps of 95% and 99%. The divergence between them is around 2%. There is no taxonomic distinction between the 13 16S sequences based on the phylogram.

One question is whether this is full length 16S? If it is there appears to be a species designation confusion in the lower clade supporting 10001 and 10013 because it is identical to to Dermococcus species and there is no way to assess which one is valid, nor whether they are really members of the same species. The very small divergence suggesting in D. abyssis might also call this into question. 10013 is in fact a robust outgroup to 10001 and the three Dermococcus spp. but the genetic distance involved is small.

However, this is not the major result. You have recovered what appears to be a large population (if that is how you sampled it) with close genetic homology to Dermococcus strain DSM20448. Given the complexity of species lineation within this genus it is impossible to say whether this represents the same species as DSM20448 but you could call it a "species complex".

Overall the population structure (if that is how it was sampled) is heavily skewed towards a bacterial species which is hitherto received minimal attention. This population skew (if that is how it was sampled), could be reversed in different ecotopes/niches/ecoregions/patients?

In summary, I wouldn't get caught up in systematics because in this genus its a bit messy. I would focus on the population structure and these results are very clear.


Standard Cutoff for Moderated T-statistics I'm looking at some microarray data. For the first time I've calculated a moderated T statistic from limma.

Is there any standard practice for where to cut off that value? For log2 fold change I usually cut off at +/- 1.5 and adjust accordingly; for adjusted p value I cut off at p<0.01 and adjust accordingly.

These moderated T statistics seem a little higher than log2 fold change and I'm not sure what is generally accepted as significant.

Thanks!

You're misinterpreting the moderated T-statistic, it's basically the fold-change divided by its variance. The p-value comes directly from that, so if you filter by moderated fold-change you're just setting an unknown (unless you go through the trouble of figuring it out) p-value threshold. Instead, do exactly as you've been doing and set appropriate (adjusted) p-value and fold-change thresholds.


How can I ask snakemake to produce a dag where each node represents a rule? I am using snakemake to create workflows. It is very convenient to visualise my DAG using snakemake --dag target{sampleA,sampleB,sampleC}File | dot -Tpdf > dag.pdf. The resulting pdf shows all the rules' dependencies to get to the target files. However, there is one node per rule/sample combination in the pdf image. When I have many samples, this creates an unreadable image with little added value. Is there an option to merge all the nodes from the same rule together so that I have the workflow for just one sample?

I think you are looking for --rulegraph instead of --dag

  --rulegraph           Do not execute anything and print the dependency graph
                        of rules in the dot language. This will be less
                        crowded than above DAG of jobs, but also show less
                        information. Note that each rule is displayed once,
                        hence the displayed graph will be cyclic if a rule
                        appears in several steps of the workflow. Use this if
                        above option leads to a DAG that is too large.


how to set database other than nr for remote blast+ search I am attempting to run a BLAST search remotely using BLAST+.

I can get search to work correctly at the command line with the following commands:

blastp -query proteins.fasta -remote -db nr -out proteins_nr.txt -outfmt 6 -evalue 1e-30

However, I would like to us the remote database titled "Microbial proteins from nr" which is what would be used for the microbial blast search accessing BLAST from the NCBI site.

I'm struggling to find the correct code to access this database instead of just "nr".

Is there a list with the codes of different databases?

Thanks in advance

I believe you're looking for env_nr? It's listed as such, under Metagenomic proteins in the blastp webpage. It appears that the word within brackets should be supplied alongside the -db parameter. A quick test with a dummy amino acid fasta file does turn up a result to a valid NCBI protein accession.

Edit: I followed up on this with a little more digging around, and it appears the database you're looking for is called Microbial_proteins. I ended up quote-searching Microbial proteins from nr on Google and exactly 2 hits turn up. The first hit, I think, is what the OP posted a screenshot of in a comment to this answer. The second is a list of databases along with keywords from a Spain-based bioinformatics firm: http://data.biobam.com/ncbi_blast_dbs_protein.pdf. I retried the dummy example with the Microbial_proteins database, and it seems to work. Again, I am not sure why I can't easily find this on an NCBI resource. Also, I am curious as to how different this database is from env_nr.