text:
Annotating conserved protein domains in bacterial assemblies
I understand the basic principles under which RPS-BLAST operates: search a protein (rpsblast) or nucleotide (rpstblastn) query against a database of protein profiles, a "reverse PSI-BLAST" search. My question has to do with the practical application of this tool to annotate conserved protein domains in bacterial assemblies.
Let's assume I have downloaded NCBI's Conserved Domain Database (CDD) and would like to annotate conserved protein domains in a bacterial genome I have just assembled from Illumina reads. I see two potential approaches.

Run a gene prediction tool on the assembly, translate the predicted gene models, use rpsblast to search the translated sequences against the CDD profiles, and then transform the matches from local protein coordinates to global genome coordinates.
Use rpstblastn to search the genome sequences against the profiles directly.

The latter option would be more compute intensive, since it requires computing a 6-frame translation of each genome sequence. But perhaps it is more sensitive in the presence of small misassemblies?
In practice, is either of these approaches preferred over the other? Are there any considerations I'm missing?

It feels like the arguments would be the same as doing blastx vs. gene prediction + blastp, and i agree with your assessment. perhaps another consideration is do you want to use this rpsblast (using PSSMs/PWMs) or something more complete like profile HMMs

The first option. Using a pre-existing tool, such as PROKKA etc.
Diversity. Bacterial diversity is staggering and species such as E. coli and Salmonella, which are closely related diverged when mammals and dinosaurs/birds diverged, so searching nucleotide sequences is a bad idea as you will get nothing. Especially as you are hoping to find distant sequences that may be missed by the commonly used annotation tools.
Operons. Bacteria, especially thermophiles, often group their protein functionally in operons. Therefore, you can use the principle of guilt by association to guess what a mysterious protein does. It is extremely powerful. Say you have Voet & Voet handy, a pathway with a hole and a mysterious protein you can guess easily what it does.
Note, if you have a protein that you are particularly interested in, I would suggest making your own PSI PSSM file. The automated system cannot read the literature or do some quick checks to see when to stop, so may be overly conservative with the clusters.


text:
Using `install.packages` with conda-managed R
I have R installed and managed by conda (miniconda) on my MacBook Pro. The version of R I use most frequently (3.5.1) is installed on the base environment and I have other version-specific environments as well. I did this as it was easier than managing R on homebrew.
The downside I now face is that I cannot use install.packages(), BiocManager::install() or biocLite() to install most packages. devtools::install_github() seems to work fine for the few packages I have installed from GitHub.
I have to use conda install -c [conda-forge|bioconda] [r|bioconductor]-package_name to install my packages and it breaks my workflow and I need to step outside of R to install packages.
Am I doing anything wrong, or is this the only way packages can be managed when R in installed through conda?
A few sample packages that I needed to install using conda:

tidyverse
tximport
deseq2
apeglm
maftools
sva

A colleague pointed me to this post: https://stackoverflow.com/q/26185978/1394178, and I tried the top 2 solutions. The Mojave solution (that installs a bunch of headers) worked for a package (refGenome) where the file wchar.h was the problem. Not sure if this would work for every package that runs into a compiler error.

What error message do you get when you try to install those bioconductor packages from R?

It's mostly C/C++ make/compile error messages. It happened so often that I stopped looking at the exact error message. I'll create a new conda env and post exact error messages for these packages.

I tried just with tidyverse. Here are the conda and R steps. The error messages in the R step are not complete, as there was only so much scroll-back history I could get to. Probably should have redirected STDERR to a file.

Installing packages outside of an R session is arguably the point of using conda.  Package management isn't something that should be done within scripts, only the loading of packages. Doing so ensures that your environment setup is somewhat robust, which should be set before you perform actual calculations.
Using conda install ... is the best way to ensure environment stability here, and, thankfully, lots of people at Bioconductor have worked very hard to make sure bioconda packages fit with the conda framework.
If you need to find packages, you can search online at https://anaconda.org, or from the command line via conda search <pkg>.

I'm really grateful to everyone creating conda install-able packages for R and bioconductor packages. I'm not looking to install packages in R scripts, only in interactive R sessions. It breaks the flow to step out of an R session and install a package. I should also mention that I do break environment stability with my habits, where my default R library is user-level (~/.Rlibs/), not environment-level. .libPaths()[2] is where the conda R library path sits.

Hi, as a follow-up: What would be the recommended way to install packages that are not available in the conda channels? Which of the other options (install.packages() or BiocManager::install() or devtools::install_github()) should be the backup.. Thanks!

Conda does handle PyPI packages well, but unfortunately it doesn't handle R packages from CRAN as nicely. Bioconductor packages available through BiocManager::install should all be listed with the name bioconductor-<packagename> in the bioconda channel, but other R packages on CRAN might not be. Personally, I've found that the easiest way to deal with CRAN-only packages is to make my own conda package using conda skeleton cran <packagename>, then install your local packge, or upload it to Anaconda cloud first before installing

source: Stack Exchange

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

sort: votes


text:
Pandas automatically rounds GWAS P-value
I am working with a specific GWAS.  If I were to run this on the command line
grep <rs_id_of_interest> GWAS.txt

I would see the GWAS p-value to be on the order of 10^-200. Very tiny p-value.
However, if I were to do the same but in pandas, with the following code,
gwas = pd.read_csv(GWAS_file,sep='\t')
gwas_rs_id = gwas.loc[gwas['marker'] == 'rs_id'].copy()

and then print gwas_rs_id, it shows the gwas p-value to be 0.0. Is there a way to read the file and not have it automatically round the tiny tiny p-values?
Thanks!

You need convert the dataframe column to a float, I'll post the code if you get stuck

I tried setting the column to np.float64, but it did not work.

The solution you likely want is here,
pd.set_eng_float_format(accuracy=x, use_eng_prefix=True)

x = whatever is required
The function set_eng_float_format has been moved around a bit and is now a top level function

You might be dealing with maximum likelihood, so convert to log likelihood -200 that sort of thing
Best idea
gwas_rs_id['logged'] = np.log(gwas_rs_id.marker)

First idea,
import decimal as D
....
gwas_rs_id['marker'] = gwas_rs_id['marker'].astype(D.decimal)

Last attempt,
pd.set_option('display.precision', x) # x whatever dp you need, this might work in combination with a log likelihood

Try a gwas_rs_id.dtype() check whether you're using float32, float64 or still in object. It should convert to scientific notion (E values) automatically.

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


text:
ChrQ definition
What is ChrQ?
I assume it then means chromosme Q or something like that yet when I google for it I get very weird answers, but as I saw that a couple times when looking on freebayes pipelines manuals I wanted to be sure.

Please [edit] your question and include an example (or a few) of where you've seen this. It's almost certainly just a place holder, as Devon suggests, but we could give a more definitive answer if you could link us to the places you've seen this used.

Answer from @devon-ryan and @terdon, converted from comments:
Are you sure it's not just a made up chromosome name that's useful in examples since people won't blindly copy/paste the command and have it work?
It's almost certainly just a place holder, as Devon suggests, but we could give a more definitive answer if you could link us to the places you've seen this used. If you only see it in relation to freebayes then yes the author is using it as “foo”.
[please update / edit to improve this answer if you can]


text:
DoHeatmap Hierarchical Clustering Seurat
I have a list of genes that I'd like to visualize using the DoHeatmap function in Seurat. However, the output of the heatmap does not result in hierarchical clustering and therefore makes it very difficult to interpret. My code is,
 DoHeatmap(object = obj, genes.use = genes),
          slim.col.label = TRUE, remove.key = T,group.label.rot = F, use.scaled = T)

Is there a way to adjust the DoHeatmap command to rank the cells by the intensity of gene expression? Here's an example output:

You need to order the marker matrix (e.g. by avg_logFC) before calling DoHeatMap.
library(dplyr)
all.markers <- FindAllMarkers(object = obj)
top20 <- all.markers %>% group_by(cluster) %>% top_n(20, avg_logFC)
DoHeatmap(object = obj, genes.use = top20$gene, slim.col.label = TRUE, remove.key = TRUE)


text:
How to interprete PlasmidSeeker results concerning plasmids?
I've been interested in plasmid detection in bacterial genomes. I tried plasmidfinder and this software found only one plasmid in each of my two bacterial genomes. I was dubious about these results, because I read many genome announcements which found something like from 6 to 10 plasmids per genome around 5 -6 MB. I then used PlasmidSeeker and what I received in the result file was a list of 9 or 11 so called "plasmid clusters", each containing a few plasmids with complete ID. I know that my genome cannot contain tens of plasmids, so my question is which one(s) of them is (are) the real one(s)?
By the way, can you suggest other bioinformatics tools specially designed to detect plasmids in bacterial genomes?


text:
How can I locate duplicated regions in a sequence?
I am facing an issue when trying to align short reads against a region in human chr5. The two Sensory Motor Neuron genes, (SMN1 and SMN2) are almost 100% identical and this causes the aligner to fail to align reads correctly since each read matches perfectly to two separate regions of the reference, leading to very low alignment scores.
If I BLAT the SMN1 gene's full genomic sequence against the human genome (hg19) on UCSC, I get (only showing the 1st two lines):
QUERY       SCORE START   END QSIZE IDENTITY  CHROM  STRAND  START       END   SPAN
--------------------------------------------------------------------------------------------------------------
NG_008691.1 35072     1 35072 35072   100.0%  chr5   +    70215768  70250839  35072
NG_008691.1 34279     1 35072 35072    99.8%  chr5   +    69340350  69375383  35034

And on the browser:

To get around this, I want to mask the SMN2 gene. However, the region of similarity extends beyond the annotated boundaries of the gene, so I would like to identify the full extent of this duplication.
So, given a large nucleotide sequence in fasta format, how can I compare it to itself and identify identical or near identical regions?

If you know approximately where are the duplicated regions, you can cut them out from the genome and align the two block using tools like MUMmer or Mauve. Although Mauve is not scalable solution (You can not scan genomes for duplicates with it) it's very handy for aligning small(ish) blocks. It can also align multiple sequences if you would like to check some tri/multiplications.
I downloaded the sequence of ch5 you linked and extracted big regions around the genes (300bp blocks). The first block is from 70100000 - 70400000 (with SMN1 sequence at 115768 - 150839) and the second from 69200000 - 69500000 (with SMN2 sequence at 140350 - 175383).
I used progressive mauve to align the two (click on File -> align with ProgressiveMauve -> add Sequence ... -> Align).

I found that the duplicated region is almost 200k long and maps between 44690 - 225520 of the first block and 19760 - 200940 of the second block, which we still need translate to the original genome coordinates, but I believe I can leave this to you :-).
P.S. you can also somehow upload annotations to Mauve, but I was lazy, so I "hand drawn" the arrow where I knew the genes are.

@KamilSJaron I can't get this to work. It keeps running out of memory even on my laptop which has 32G, of which 20 were available. I even tried aligning 2 files of only 100bp each and it crashed. I managed to get it to run on a server with 41G of RAM available, but I can't get the graphical output that way. I see you ran it on macOS, so it seems like a bug of the Linux version :(

GeneBlocks is a Python library for comparing DNA sequences. It can be used to find common blocks (identical regions) in a group of DNA sequences. Using a modified version of your example (I duplicated a small region):
import geneblocks

sequence = geneblocks.load_record('sequence.gb', name='seq')
sequence2 = geneblocks.load_record('sequence.gb', name='seq2')
blocks = geneblocks.CommonBlocks.from_sequences([sequence, sequence2])

blocks.plot_common_blocks()

Note that it wasn't designed to plot only one sequence, so as a workaround the sequence is given twice, with different names. As an alternative, you can use one sequence and plot with DNA Features Viewer:
import dna_features_viewer

blocks = geneblocks.CommonBlocks.from_sequences([sequence])
seq_to_plot = blocks.sequences_with_annotated_blocks()["seq"]
graphic_record = dna_features_viewer.BiopythonTranslator().translate_record(seq_to_plot)
ax, _ = graphic_record.plot(figure_width=10, strand_in_label_threshold=7)

See the documentation for more details. A web interface is also provided at EGF CUBA: find common blocks.
Disclaimer: I'm the current maintainer of GeneBlocks


text:
How to select only RNA with Hetero atoms from pdb file with python?
I'm trying to separate RNA from protein in a complex protein/RNA PDB file and I want all RNA info with the hetero atoms in between the bases BUT without H20 etc. In short I want RNA part of pdb file without discontinuous lines.
I managed to separate RNA from protein with Bio PDB Select but it consider hetero atoms as amino acid when I use is_aa(residue). So hetero atoms wont appear in my "only RNA" file.
from Bio.PDB import *
from Bio.PDB import PDBParser, PDBIO, Select
import os

class ProtSelect(Select):
    def accept_residue(self, residue):
        return 1 if is_aa(residue) == True else 0

class RNASelect(Select):
    def accept_residue(self, residue):
        return 1 if is_aa(residue) == False and residue.id[0] != "W" else 0

pdb = PDBParser().get_structure("2bh2", "pdb2bh2.ent")
io = PDBIO()
io.set_structure(pdb)
io.save("seqprotest.pdb", ProtSelect())
io.save("seqRNAtest.pdb", RNASelect())

EDIT: I managed to do it:
class RNASelect(Select):
    def accept_residue(self, residue):
        return 1 if ((residue.id[0] != 'H_ MG') and ((residue.get_resname() == '  A' or residue.get_resname() == '  U' or   residue.get_resname() == '  G' or  residue.get_resname() == '  C' or residue.get_resname() == '  T') or  (residue.id[0][0:2] == 'H_' and (residue.get_resname()[2] == 'U'  or residue.get_resname()[2] == 'A' or residue.get_resname()[2] == 'C' or residue.get_resname()[2] == 'G' or residue.get_resname()[2] == 'T'))))  else 0

It's not pretty but it do the job.
The idea is to get everything that finish by a nucleotide. I figured out that hetatom's names of RNA chains always finish by a nucleotide letter so i just have to check everything that begins by 'H_' and ends by 'A,T,C,G,U'.
Just remove 'H_ MG' and you're good.

I've not done , RNA vs protein but separating protein vs protein is easy in Pymol. Whether Pymol will separate heteroatoms I have no idea

I don't remember now but when I had to work with PDB the biopython tutorial had a long section about it. I think there was some other property/method to find if it was really an amino acid or something else. It is worth exploring

@Raph post the solution as an answer to help other too (It can be missed if it is in the body of the question)

Solution I found: (c/p of body)
class RNASelect(Select):
    def accept_residue(self, residue):
        return 1 if ((residue.id[0] != 'H_ MG') and ((residue.get_resname() == '  A' or residue.get_resname() == '  U' or   residue.get_resname() == '  G' or  residue.get_resname() == '  C' or residue.get_resname() == '  T') or  (residue.id[0][0:2] == 'H_' and (residue.get_resname()[2] == 'U'  or residue.get_resname()[2] == 'A' or residue.get_resname()[2] == 'C' or residue.get_resname()[2] == 'G' or residue.get_resname()[2] == 'T'))))  else 0

It's not pretty but it do the job. The idea is to get everything that finish by a nucleotide. I figured out that hetatom's names of RNA chains always finish by a nucleotide letter so i just have to check everything that begins by 'H_' and ends by 'A,T,C,G,U'. Just remove 'H_ MG' and you're good.


text:
for loop cutadapt for files on a single directory
I have a bigfolder where i have a lot of fastq.gz files and I want to remove the adapters from all of them.
I am trying then the following loop:
for file in LotsofRaws\/cannam/*_R1.fastq.gz; do cutadapt
     -a runs/new/Adapters/IlluNexteraPE.fa
        -o LotsofRaws/RemovedAdapter/Adapted$file
           $file done;  done

However I am getting the following error:
cutadapt: error:

*It appears you want to trim paired-end data because you provided two input files, but then you also need to provide two output files (with -o and -p) or use the --interleaved option.*

But I just want to take the adapters from the files I am not looking to do it in the same way I would do with a single file is there any way to do this because if I try to do --interleaved I get:
cutadapt: error: Character '/' in adapter sequence 'RTNS/NEW/ADAPTERS/ILLTNEXTERAPE.FA' is not a valid IUPAC code. Use only characters XACGTURYSWKMBDHVN.

This is cutadapt 2.3 with Python 3.6.7
and if instead of using the file I directly use the sequences:
for file in LotsofRaws/cannam/*;
         do cutadapt -a CTGTCTCTTATACACATCT
               -a AGATCGGAAGAGC
               --interleaved
               -o LotsofRaws/RemovedAdapter/$file   $file ; done

I get:
FileNotFoundError: [Errno 2] No such file or directory: 'LotsofRaws/RemovedAdapter/LotsofRaws/cannam/NAM7_Glenlea_R1.fastq.gz'
During handling of the above exception, another exception occurred:

many times
What should I then do for general variant calling it is appropriate to do the --interleaved option?
Should I just the do the trimming manually or can I loop it in any way that still resembles doing it one by one but without me having to go and write line and line again.
So the files on the folder look like:
NAM11_BW415_R1.fastq.gz         NAM13_BW968_R2.fastq.gz     NAM17_AC_Barrie_R1.fastq.gz  NAM3_PI648600_R2.fastq.gz  NAM7_Glenlea_R1.fastq.gz   PI648598_R2.fastq.gz
NAM11_BW415_R2.fastq.gz         NAM13_Frontana_R1.fastq.gz  NAM17_AC_Barrie_R2.fastq.gz  NAM3_PI648842_R1.fastq.gz  NAM7_Glenlea_R2.fastq.gz   PI648611_R1.fastq.gz
NAM11_Reeder_R1.fastq.gz        NAM13_Frontana_R2.fastq.gz  NAM17_BW388_R1.fastq.gz      NAM3_PI648842_R2.fastq.gz  NAM7_Laura_R1.fastq.gz     PI648611_R2.fastq.gz
NAM11_Reeder_R2.fastq.gz        NAM14_BW874_R1.fastq.gz     NAM17_BW388_R2.fastq.gz      NAM4_PI613279_R1.fastq.gz  NAM7_Laura_R2.fastq.gz     PI648672_R1.fastq.gz
NAM11_Webill_R1.fastq.gz        NAM14_BW874_R2.fastq.gz     NAM17_RedFife_R1.fastq.gz    NAM4_PI613279_R2.fastq.gz  NAM7_Marquis_R1.fastq.gz   PI648672_R2.fastq.gz
NAM11_Webill_R2.fastq.gz        NAM14_Glenn_R1.fastq.gz     NAM17_RedFife_R2.fastq.gz    NAM4_PI613280_R1.fastq.gz  NAM7_Marquis_R2.fastq.gz   PI648724_R1.fastq.gz
NAM12_BW431_L002_R1.fastq.gz    NAM14_Glenn_R2.fastq.gz     NAM1_PI648517_R1.fastq.gz    NAM4_PI613280_R2.fastq.gz  NAM8_FL62R1_R1.fastq.gz    PI648724_R2.fastq.gz
NAM12_BW431_L002_R2.fastq.gz    NAM14_Unity_R1.fastq.gz     NAM1_PI648517_R2.fastq.gz    NAM4_PI648524_R1.fastq.gz  NAM8_FL62R1_R2.fastq.gz    PI648768_R1.fastq.gz
NAM12_BW431_L008_R1.fastq.gz    NAM14_Unity_R2.fastq.gz     NAM1_PI648519_R1.fastq.gz    NAM4_PI648524_R2.fastq.gz  NAM8_Muchmore_R1.fastq.gz  PI648768_R2.fastq.gz
NAM12_BW431_L008_R2.fastq.gz    NAM15_BW881_R1.fastq.gz     NAM1_PI648519_R2.fastq.gz    NAM5_PI648523_R1.fastq.gz  NAM8_Muchmore_R2.fastq.gz  PI648860_R1.fastq.gz
NAM12_Penhold_L002_R1.fastq.gz  NAM15_BW881_R2.fastq.gz     NAM1_PI648541_R1.fastq.gz    NAM5_PI648523_R2.fastq.gz  NAM8_Stettler_R1.fastq.gz  PI648860_R2.fastq.gz
NAM12_Penhold_L002_R2.fastq.gz  NAM15_Lillian_R1.fastq.gz   NAM1_PI648541_R2.fastq.gz    NAM5_PI648585_R1.fastq.gz  NAM8_Stettler_R2.fastq.gz  PI648862_R1.fastq.gz
NAM12_Penhold_L008_R1.fastq.gz  NAM15_Lillian_R2.fastq.gz   NAM2_PI648493_R1.fastq.gz    NAM5_PI648585_R2.fastq.gz  PI277012_R1.fastq.gz       PI648862_R2.fastq.gz
NAM12_Penhold_L008_R2.fastq.gz  NAM15_PT245_R1.fastq.gz     NAM2_PI648493_R2.fastq.gz    NAM5_PI648622_R1.fastq.gz  PI277012_R2.fastq.gz       PI648864_R1.fastq.gz
NAM12_Sumai3_L002_R1.fastq.gz   NAM15_PT245_R2.fastq.gz     NAM2_PI648602_R1.fastq.gz    NAM5_PI648622_R2.fastq.gz  PI648478_R1.fastq.gz       PI648864_R2.fastq.gz
NAM12_Sumai3_L002_R2.fastq.gz   NAM16_BA51B92_R1.fastq.gz   NAM2_PI648602_R2.fastq.gz    NAM6_Neepawa_R1.fastq.gz   PI648478_R2.fastq.gz       run_701_150904_7001410_0194_BC79VNANXX_NAM.list
NAM12_Sumai3_L008_R1.fastq.gz   NAM16_BA51B92_R2.fastq.gz   NAM2_PI648609_R1.fastq.gz    NAM6_Neepawa_R2.fastq.gz   PI648487_R1.fastq.gz       run_760_160912_7001343F_0065_AC9EV0ANXX_NAM.list
NAM12_Sumai3_L008_R2.fastq.gz   NAM16_BW928_R1.fastq.gz     NAM2_PI648609_R2.fastq.gz    NAM6_PI648532_R1.fastq.gz  PI648487_R2.fastq.gz
NAM13_BW462_R1.fastq.gz         NAM16_BW928_R2.fastq.gz     NAM3_PI648511_R1.fastq.gz    NAM6_PI648532_R2.fastq.gz  PI648566_R1.fastq.gz
NAM13_BW462_R2.fastq.gz         NAM16_BW963_R1.fastq.gz     NAM3_PI648511_R2.fastq.gz    NAM6_Thatcher_R1.fastq.gz  PI648566_R2.fastq.gz
NAM13_BW968_R1.fastq.gz         NAM16_BW963_R2.fastq.gz     NAM3_PI648600_R1.fastq.gz    NAM6_Thatcher_R2.fastq.gz  PI648598_R1.fastq.gz

I see a number of issues with your code: 1. Why does the first loop have a \  after LotsOfRaws? 2. In your second loop, you're passing in the entire file path as $file, whereas you should be setting something like file_bn=$(basename $file) and using $file_bn instead of $file with the -o option.

You have various different issues here. First, the first version of the command you show would give a syntax error, not the error you show, so I am guessing what you actually ran was this:
for file in LotsofRaws\/cannam/*_R1.fastq.gz; do
    cutadapt -a runs/new/Adapters/IlluNexteraPE.fa \
             -o LotsofRaws/RemovedAdapter/Adapted$file $file done;
done

The error there is that you have the word done after $file. So what cutadapt sees is:
cutadapt [OPTIONS] $file done

So it reads that as being given two input files instead of one. What you wanted to do was this:
for file in LotsofRaws\/cannam/*_R1.fastq.gz; do
    cutadapt -a runs/new/Adapters/IlluNexteraPE.fa \
             -o LotsofRaws/RemovedAdapter/Adapted$file $file
done

Or, on one line, this:
for file in LotsofRaws\/cannam/*_R1.fastq.gz; do cutadapt -a runs/new/Adapters/IlluNexteraPE.fa -o LotsofRaws/RemovedAdapter/Adapted$file $file; done

One done, not two and separated from the command by a newline or a ;.
Your next attempt, with --interleaved failed because the -a option expects an actual sequence, not a file unless you use it with file:. See cutadapt --help (emphasis mine):

Finding adapters:
  Parameters -a, -g, -b specify adapters to be
  removed from each read (or   from the first read in a pair if data is
  paired). If specified multiple   times, only the best matching adapter
  is trimmed (but see the --times   option). When the special notation
  'file:FILE' is used, adapter sequences   are read from the given FASTA
  file.

So, the -a expects to be given a sequence unless you use it as -a file:/path/to/file. Because you used it as -a runs/new/Adapters/IlluNexteraPE.fa instead of -a file:withruns/new/Adapters/IlluNexteraPE.fa, it tried to interpret runs/new/Adapters/IlluNexteraPE.fa as a sequence and complained because the character / isn't valid in a DNA sequence.
Finally, your last lop failed because you're giving the wrong name. The first rule of debugging is "print all the things". When you see a command that doesn't work as you expect, your first step should always be to print out the command and see what is actually being run. When you do a for file in /long/path/to/file loop, the entire path is included in your $file variable. Try running this:
for file in LotsofRaws/cannam/*; do
    echo "cutadapt -a CTGTCTCTTATACACATCT -a AGATCGGAAGAGC --interleaved -o LotsofRaws/RemovedAdapter/$file $file";
done

You will see that that will print commands like this:
cutadapt -a CTGTCTCTTATACACATCT -a AGATCGGAAGAGC --interleaved \
         -o LotsofRaws/RemovedAdapter/LotsofRaws/cannam/NAM11_BW415_R1.fastq.gz \
          LotsofRaws/cannam/NAM11_BW415_R1.fastq.gz

What you need is to take the name of the file, without the path. And you don't want to escape the slash after LotsofRaws. So putting all this together and going back to your first attempt, what you want is something like this:
for f in LotsofRaws\/cannam/*_R1.fastq.gz; do
    fname=$(basename file)
    cutadapt -a "file:runs/new/Adapters/IlluNexteraPE.fa"
             -o "LotsofRaws/RemovedAdapter/Adapted${fname}" $file
done


text:
What are the right parameters to trim a small RNA transcriptome with trimmomatic?
I'm having some problems with finding the right parameters to trim my small RNA Illumina reads (51 nt long) with Trimmomatic.
Before trimming, one of the samples (21M reads) looks like this:

So for my understanding, it's quite good and the last two issues can be easily solved with adapter trimming.
I then tried to trim it with the command
java -jar /opt/Trimmomatic-0.38/trimmomatic-0.38.jar SE -threads 24 -phred33 FemaleMito3.fastq FemaleMito3_trimmed.fastq ILLUMINACLIP:adapters.ultimate.fa:2:30:10 AVGQUAL:25 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:18

but I got the following result:

So now I have a problem with the sequence quality per tile (why?), an overall drop of the quality at the end of the reads and a warning with the sequence length distribution.
Then I tried trimming just the adapters with this command:
java -jar /opt/Trimmomatic-0.38/trimmomatic-0.38.jar SE -threads 24 -phred33 FemaleMito1.fastq FemaleMito1_noadapters.fastq ILLUMINACLIP:adapters.ultimate.fa:2:30:10 AVGQUAL:25

but again I get some new issues:

Again the per tile sequence quality issue, and now a warning about the GC content.
I didn't find anywhere some precise instruction on how to trim a small RNA
library with Trimmomatic. How should I properly trim my small RNA library with Trimmomatic?

Relevant: https://bioinformatics.stackexchange.com/a/6788/96.

The command where you trim with adapters and by quality is perfectly fine. That FastQC isn't perfectly happy is expected. It's a tool made for whole genome sequencing QC, you SHOULD see a number of "fail"s with any RNA-seq protocol. Steps that should always fail in RNA-seq:

per-base sequence content (there's "random priming" performed that isn't completely random, so the first ~8 bases show a typical bias pattern).
per-sequence GC content (particularly in sRNA-seq you're enriching for a subset of genes, they won't have some nice GC distribution and that's perfectly OK).
Sequence duplication level (if this doesn't fail then everything is uniformly expressed...which would be quite unlikely)
Overrepresented sequences (as above, though this will be worse for sRNA-seq since reads will often cover the entire transcript, so there's even less sequence diversity)

Thank you for your answer. I expected FastQC not being happy with the sequence duplication level and overrepresented sequences, but I couldn't explain the other two. What about the "per tile sequence quality"? Did adapters inflated artificially the quality of those tiles? Should I do something about that (filter those reads out...) or can I leave them as is?

Yes, it appears that the adapters inflated quality a bit. Don't worry about failed clusters, there's no point in trimming that minor an issue off.


text:
Proteomics data Vs Transcriptomics data?
I want to use either of Proteomics or Transcriptomics data for integrating it into my kinetic model. Before proceeding, I want to know what are the advantages of using either of them so that I could make an informed decision on it!
Many studies have shown that the best we can do is integrating both transcriptomics and proteomics data with our kinetic model, but I've some time constraints and have to proceed with only one of those.
My effort and findings: I've found from discussions with researchers that gathering transcriptomics data has an amplification step which increases the chance of finding a particular one whereas gathering proteomics data has no such step but has fragmentation and then rejoining which creates many problems(such as splice variants etc) and thus leads to a loss of data. But a PostDoc told me that even after the loss of that data, I'll get more information from Proteomics data.
I want to know such type of points and want to know if these are valid or not!

You want to compare two different techniques but do you have data from both of them for your model? Theoretically they can be whatever, but then in your experiment you can get unusable data. Also it seems like you have two implementations of the model, so it might be influence the accuracy of the model. Did you try the model with simulated data of both types? (or with freely available information on the web)

I am yet to gather the data(either transcriptomics or proteomics) for my model. I have the list of genes ready for which I want it. I have only one implementation of the model only. I haven't tried the simulation with any data as I don't have it yet. I want to make an informed decision on it. Please suggest @llrs

I did some research and have put down my findings. Please add any suggestions that you want. @llrs

Ultimately, proteomics data can actually measure the actual concentrations of the enzymes that you are interested, while transcriptomics only measures their rates of change due to ongoing transcription from mRNA.
Rates of change are often treated as proxy for concentrations. This makes sense when you have a constant degradation/dilution environment and are near steady state, since then the rate of change is proportional to the converged concentration. If you're interested in more dynamic environments or believe that modulation of degradation or sequestration might play a significant role in a system, however, that assumption breaks done.
Thus, if you have to choose between proteomics and transcriptomics data, I would suggest proteomics.
As a final note, however, since the link suggests you're interested in metabolites, however, have you consider going directly to metabolomics data?  That would let you really get direct information about the particular species of interests.

This is what I found from doing some research. Comments are welcome at any point!

Capture percentage for data gathering:

Transcriptomics Data: There's an amplification step in Transcriptomics data gathering methods. Hence, it's possible to capture almost the totality of the Transcriptome using those methods(scRNA seq methods, Nanopore tech, Spatial transcriptomics, etc)
Proteomics Data: There's no amplification step in Proteomics data gathering methods. Those methods have fragmentation and a rejoining step which leads to a loss of the data.

Scalability of methods:

Transcriptomics Data: Scalable methods available
Proteomics Data: Less scalable methods for protein studies

Reference availability:

Transcriptomics Data: Fully annotated references available on consortiums such as ensembl biomart, etc
Proteomics Data: Universal and Comprehensive Human Proteome reference is still in question and because of the incomplete and inaccurate references, much of the newly generated data is being rendered useless

Uniformity

Transcriptomics Data: Lots of publications on pipelines, and their benchmarking available
Proteomics Data: Lack of uniformity across labs/research groups and lack of related literature lead to differences in protein fragmentation and solubilization and differences in algorithms to run analyses

Technical Bias

Transcriptomics Data: Many methods available to tackle cell-bias, noise in data, and batch effects, so as to capture maximum biological variability
Proteomics Data: Cell-wide Mass spectrometry has a bias towards the identification of peptides with higher concentration or contamination from other experiments

On the plus side, proteomics is actually measuring what's most relevant to your model.

@DevonRyan yes, sometimes I think that the field is focused too much on transcriptomics... I hope to see more usage (and learn more) on proteomics. On which side do you lean Akhil? BTW Did you take into account the experimental and bioinformatic experience as well as the cost of each to obtain (and store/process) the data you need? (just throwing ideas that are worth considering from a practical point of view)

I'm more inclined to use the proteomics data only but my PI has told me to go with whatever data will be of higher quality as then we can draw some inferences more confidently. Also, I didn't get your last point. Can you please elaborate on that? @llrs

Is anyone on your work (lab?) that has done RNA extractions for RNAseq or anyone that have done the preparation for mass spectrometry? Assuming you have them. How much would it cost to have the same quality of data from both techniques? Once you have the data, is there anyone that has the knowledge and resources to process it? In RNAseq you need to align to the reference (computationally expensive and time consuming), normalize/prepare it for your model while store all the data and the original files you got the fastq.gz. How do you plan to do that?

My lab is purely a computational one. We don't have the resources or people to perform the said experiments. So, my PI has asked me to gather the data from the existing literature. In the same way, we gathered the metabolites and flux data from the existing literature on Mass spectrometry and have used it to construct a kinetic model. Just to tell you I'm very well familiarised with the reference alignment, normalization, etc as I've done an internship in the scRNA seq analysis, but I don't think it'll be useful to me here in anyway. @llrs


text:
How should I deal with segmental duplications when aligning NGS reads to a reference genome?
This is a follow-up of my other question. I have been having trouble calling variants in the human SMN1 and SMN2 genes, because the human genome has a large segmental duplication there and these two genes (and their surrounding sequences) are essentially identical :

Click here to see the original image on the UCSC browser.
Since variants in the SMN1 gene have been linked to disease, it is important for me to be able to call them in my NGS data. However, because of this segmental duplication, aligning reads to this region correctly is very hard since the presence of the duplication makes each read align almost perfectly to separate locations in the genome, leading to low alignment scores. The low scores will, in turn, block variant callers from calling in those regions.
I thought that the only real solution for this would be to mask the SMN2 gene and its associated duplication in my reference genome which would allow the aligner to align all reads to the region of SMN1. I would then miss any real variants in SMN2, but since the primary gene of interest for  spinal muscular atrophy is SMN1 and not SMN2, I would be OK with that.
My questions are:

Is this a reasonable approach? I was told that "generally, masking segdups is not recommended" by someone who most certainly knows what he's talking about. So is this not a good approach for my stated goal?

If not, then what else can I do that will let me call variants in the SMN1 gene?

I would like something similar to the approach described here: mask one of the regions to rescue mutations in the other.

Short answer for short reads: you can't. Not reliably, at least.
There's a known error mode with aligning short reads to reference genomes where it's quite good at mis-calling reads within repetitive regions. See, for example, this paper, which compares Illumina vs Nanopore for detecting bacterial strains.
If you know that a region is likely to contain a segmental duplication and/or deletion (e.g. due to unexpected changes in read coverage), then it'd be a good idea to double-check any discovered variants (or unexpected non-variants) using a different technology that does not have the same coverage issues. I recommend a long-read sequencing technology that can handle regions at least as large as the segmental duplication of interest (e.g. PacBio up to ~15kb, Nanopore for anything longer), so that you have a high chance of mapping reads uniquely.

I agree that a longer read technique would help but unfortunately that isn't an option for me as I am processing the data generated by others. I need to improve the results of using short reads and don't have the option of using long. But I don't see why you say "you can't". It turns out I most certainly can, and the approach described in the paper you linked to in the question (thanks, I should have included that) works. The complication is dealing with different repeat numbers and ploidy, but the general approach certainly seems to work so far and is giving more accurate results.


text:
Which statistical test to use for comparing two groups?
I have data for 396 transcriptional factors (TF), where I have calculated the number of Transcriptional targets for each transcriptional factor present between the two groups:
Group-I(G-I): where co-exp genes have a correlation coefficient is more than 0.6 and
Group-II(G-II): where the co-exp genes have r values less than -0.6.
I have co-expression map for all the human transcription factors. My hypothesis is that transcriptional targets of a transcriptional factor are among the highly co-expressed gene for that transcription factor. For that, I have divided my group into high and low co-exp genes.
My question which statistical test I should use for evaluating the difference between the two groups?
  TF    G-I G-II
  AATF  1   0
  ABL1  5   1
  AES   0   1
  AHR   1   2
  AIP   2   0
 APBB1  0   1
APEX1   2   1
  AR    0   2
ARNT    3   0

Post a histogram or similar of the two groups, it would help to see if they have a particular distribution.


text:
How to find novel transcripts using GFFcompare?
I am trying to find novel transcripts from an RNA-seq database. Based on the advice I got, it seemed that using Stringtie for transcript assembly is a good way to go, and it supports novel transcript discovery even with the reference GTF file provided (which apparently significantly improves the assembly process). So I tried to follow the protocol provided in the Nature Protocols paper that described the use of Stringtie; the paper also suggested using GFFcompare for comparing the assembled transcripts with the reference, to quantify and find the novel transcripts.
I followed the entire protocol, and this was the output of GFFcompare when comparing the reference GTF with that of one of the samples after mapping -
#= Summary for dataset: ./hisat2/ERR188044_chrX.gtf
#-----------------| Sensitivity | Precision  |
        Base level:    51.7     |    79.5    |
        Exon level:    46.7     |    85.2    |
      Intron level:    47.2     |    93.9    |
Intron chain level:    31.2     |    64.4    |
  Transcript level:    31.6     |    52.3    |
       Locus level:    36.6     |    50.1    |

     Matching intron chains:     580
       Matching transcripts:     664
              Matching loci:     397

          Missed exons:    4395/8804    ( 49.9%)
           Novel exons:     426/4874    (  8.7%)
        Missed introns:    3832/7946    ( 48.2%)
         Novel introns:      83/3992    (  2.1%)
           Missed loci:     610/1086    ( 56.2%)
            Novel loci:     273/795 ( 34.3%)

 Total union super-loci across all input datasets: 749
1270 out of 1270 consensus transcripts written in 188044_only.annotated.gtf (0 discarded as redundant)

In the paper however, the number of novel genes and transcripts is clearly mentioned -

So, I was wondering how they were able to calculate the number of novel transcripts and genes, from the data that GFFcompare provides (they even mention in the paper that they got it from GFFcompare). I see the numbers related to novel exons and novel introns, but how do I translate it to the number and identities of the novel transcripts/genes?
Cross-posted in Biostars. Let me know if I shouldn't do this.

The output of gffcompare includes several files per run (just like cuffcompare). Example for a run:
$ ls cuffcmp.* | sed 's/\t/\n/'
cuffcmp.combined.gtf
cuffcmp.loci
cuffcmp.output.gtf.refmap
cuffcmp.output.gtf.tmap
cuffcmp.stats
cuffcmp.tracking

From my experience, the easiest file to manipulate is the tmap file:

Tab delimited file lists the most closely matching reference transcript for each query transcript. There is one row per query transcript, and the columns are as follows:
1     Reference gene name    The gene_name attribute of the reference GTF record for this transcript, if present. Otherwise gene_id is used.
2     Reference transcript id     The transcript_id attribute of the reference GTF record for this transcript
3     Class code  The type of relationship between the query transcripts in column 4 and the reference transcript (as described in the Class Codes section below)
4     Query gene id   The query (e.g., Stringtie) internal gene id
5     Query transcript id     The query internal transcript id
6     Number of exons     The number of exons in the query transcript
7     FPKM    The expression of this transcript expressed in FPKM
8     TPM     the estimated TPM for the transcript, if found in the query input file
9     Coverage    The estimated average depth of read coverage across the transcript.
10    Length  The length of the transcript
11    Major isoform ID    The query ID of the gene's major isoform
12    Reference match length  The length of the longest overlap with a reference, ‘-’ if there is no such exonic overlap

The important column to look at for transcript classification is actually the 3rd, "class code". The content is meaningful if you used option -r when running gffcompare.
In the tmap file you can access info on both transcript and gene ID directly, both from the assembler and the reference annotation.
Now what you need is:

n. assembled genes
n. novel genes
n. novel transcripts
n. transcripts matching annotation

At this point, it is essential to understand what the authors mean by "novel" (gene/transcript), does it mean "non previously annotated" (= class code u, for "unknown")? Also, what about "matching", does it mean "perfectly matching" (= class code =, for "complete, exact match of intron chain")?
Under these assumptions, you can do some easy stats with some bash help.

n. assembled genes:
$cat cuffcmp.output.gtf.tmap | sed "1d" | cut -f4 | sort | uniq | wc -l
n. novel genes:
$ cat cuffcmp.output.gtf.tmap | awk '$3=="u"{print $0}' | cut -f4 | sort | uniq | wc -l
n. novel transcripts:
$ cat cuffcmp.output.gtf.tmap | awk '$3=="u"{print $0}' | cut -f5 | sort | uniq | wc -l
n. transcripts matching annotation:
$ cat cuffcmp.output.gtf.tmap | awk '$3=="="{print $0}' | cut -f5 | sort | uniq | wc -l

I hope this helps. I think it's an easier approach than trying to "translate the numbers related to novel exons and novel introns to the number and identities of the novel transcripts/genes" from the stats file.
Although you have to think a bit what "class codes" to include in each count. Maybe my assumptions are too strict for your use case.

Thank you so much for this detailed response, this was life-saving. I had been trying to find a way out for so long on multiple forums! Turns out I just need to look better at the files that the tool had already generated for me. Whats interesting is that using the code you describe, I get about ~80K assembled genes and ~220K transcripts matching annotation! I checked why this could be and found that most entries in the column for assembled gene (query gene ID) are repeated, while there are no repeats in the column for the query ID (used for transcripts matching annotation). I wonder why!


text:
Alternative to enrichR for enrichment analysis?
I have been using EdgeR from the Bioconductor package to perform differential expression analysis on bulk RNASeq data sets.
This analysis gives me an output dataframe that has a list of all the genes, ranked by the their adjP values. This dataframe also contains the t-statistic and lfc values for each gene as well.
I have recently been using enrichR to perform gene set enrichment analysis on all the differentially expressed genes to determine what biological, molecular, transcriptional, etc pathways these genes are involved in. I really like the enrichR package for a few key reasons:
1.) It has access to many different databases that I can call on with a simple vector of their names. Most of my data is murine RNASEq, so for those analyses, I can simply select the murine databases:
enrichR_dbs <- c("WikiPathways_2019_Mouse", "KEGG_2019_Mouse", "GO_Biological_Process_2018", "GO_Cellular_Component_2018", "GO_Molecular_Function_2018", "GWAS_Catalog_2019", "TRRUST_Transcription_Factors_2019")

2.) It uses a single function, enrichR, to perform the enrichment. The arguments for that function are the list of gene symbols and the list of databases.
enrichr(as.character(Gene_Dataframe$SYMBOL[1:500]), enrichR_dbs)

Anyway, I am looking for another enrichment analysis because I want to be able to see if I get similar pathways for my genes when using another method/algorithm.
Unfortunately, I can't seem to find another enrichment analysis( competitive or self-contained) that is both as user-friendly and that gives me easy access to murine databases (without requiring to convert my genes through changes in capitalization, like you need to with GSEA) as enrichR.
Does anyone have recommendations about enrichment analysis R packages that don't require too many more inputs than the dataframe I generated through my EdgeR analysis?
P.S. I'm ok with downloading databases if I need to. Though inclusion of a PDF with a tutorial would be fantastic, if possible.

A couple of years ago I wrote a blog post about GSEA methods in Bioconductor, in addition to these I would recommend clusterProfiler or ReactomePA.  In most of them you need to provide the gene set/pathway database.
If you want to do gene ontology enrichment then I recomment the topGO package.
In most, if not all, of these methods you just need an ordered list and a gene set source.

Try gProfileR. It is very intuitive and easy to use in R.


text:
rentrez to retrieve NIH files
I'm trying to use rentrez package to retrieve NIH files, and ran the following code (I followed what was on the rentrez tutorial):
pubmed_search <- entrez_search(db = "pubmed",
                    term = "Case Reports[Filter] AND cardiovascular disease AND English[lang] AND 2009:2019[PDat])",
                    retmax = 792711,
                    use_history = TRUE)

for(i in seq(1,10000,100)){
    recs <- entrez_summary(db="pubmed",
                           web_history=pubmed_search$web_history,
                           retmax=100,
                           retstart=i)   cat(
seq_start+99, "sequences downloaded\r") }

However, instead of the 10,000 texts I expected, I only downloaded 100 files. Could so

Welcome to Bioinformatics.SE! It looks like you posted the question before you were done typing. No worries, just use the edit button to update it. Also, please indent code with 4 spaces or wrap it in triple backticks (\```) so that it displays in a legible way.

Also, include an API code in your Pubmed search script to increase the query limit.

I don't program in R, but replace the last line of you code with the following and I suspect you will find it works,
for(i in seq(1,10000,1)){ recs <- entrez_summary(db="pubmed", web_history=pubmed_search$web_history, retmax=100, retstart=i) cat(seq_start+99, "sequences downloaded\r") }

Rationale I would guess seq is similar to the linspace function in Python, and put simply 10000 (total sample size) / 100 (the step size) = 100, 100 being the number of downloads you observed. In the code above its 10000/1. I would also suspect "step size = 1" is default, if this was correct then seq(1,10000) would give you your required output.
As a random aside I do like the presence of curly brackets defining the loop, it reminds me of Perl.


text:
Statistical test design question
I have condition as such mouse liver and spleen data.
So in liver I have WT (Wild Type) and KO (knock out) as well as in spleen WT (Wild Type) and KO (knock out).
I have to compare the parasite load in WT and KO condition in both liver and spleen.
In terms of data points what I got is this:

10 mice for WT, 6 mice for KO in case of liver
10 mice for WT, 5 mice for KO in case of spleen.
WT 20 8 21 9 22 10 23 11 24
KO 12 25 26 27 28

So far what I know is I have to perform a parametric test.
But I'm not sure if I have to go for a paired t-test or an unpaired t-test.
Because here I'm comparing with between two groups, only so I don't have to go for anova if I'm not wrong. Secondly there is a difference in number of observation that is available for WT as compared to KO. So would that affect my choice of test as well.
Suggestion would be highly appreciated

Please take some time to clarify your question. It's unclear what you're ultimately asking. And your description of the data is unclear as well. For example, what you mean by days? Is this temporal data (WT parasite load at 1, 2, ..., 10 days), or 10 WT replicates for (in which case the significance of "days" is unclear).

sorry will edit the question give me a moment

edited the question and given you a sample scenario

Your edits help, but the post is still very verbose and unclear. How/Why do you know you have to perform a parametric test? What do you want your test to tell you?

". How/Why do you know you have to perform a parametric test?" i guess if the data is normally distributed then parametric if not then non-parametric . "What do you want your test to tell you" I want to see there is difference in the parasitic load between WT and KO cases one in the liver and other in the spleen

Here's a combination of general advice and my thoughts about your specific data:
If the samples are not matched, e.g. by time point, then it will be unpaired. You need a way to say sample 1 KO and sample 1 WT are uniquely matched to be paired, sample 2 KO and sample WT will be paired etc ...  pairing is great for a time series. If you have replicates within the paired groups then ANOVA is the way to go, particularly if the sample sizes are uneven.
Note: in your case the sample sizes are uneven, you cannnot perform a paired T-test between two uneven sample sizes. HOWEVER you appear to a natural pairing, but I don't know whether this is biologically justified, i.e. liver and spleen. So you appear to have a biological pairing with uneven sample sizes. If this is correct it suggests a 1-way ANOVA is a better idea than a T-test.
Parametric test: check for normality in each data set, e.g.

A Q-Q plot against the normal distribution (a cool way to do this)
Homogeneity of variances is Hartley's Fmax test (good for T-test and ESSENTIAL for 1-way ANOVA)
68% of the data should be less than ±1 standard deviation around the mean, 95% less than ±2 standard deviation and 99% less than ±3
Kolmogorov-Smirnov and/or Shapiro-Wilk test
Skewness & kurtosis magnitudes < 1.96 respectively

Note if you combine liver and spleen data sets and find a nice normal distribution, and one or both is not normal by itself, then you would need to consider the biological justification for combining the data sets.
If the data passes all those tests, you are good to go. If not .... gets complicated,

Normalisation, this is where you transform the data, in your case I'm not sure, but typically log transformation and then perform the checks for normality above. There are lots of ways to normalise data.
Repeat the T-test/1-way ANOVA with normalised data
Ignore it - yup its top level statistical advice, it is context specific and never a universally shared idea ... I suspect that is where the advice "just use a parametric test" came from. Please be aware if you do this, someone could turn around and sau "you've done it wrong". The statistical arguments about why this assumption can be valid get complicated.
Use a non-parametric test, orthodox statistical advice, in this case it is the Mann-Whitney U test and Kruskal–Wallis test if you have paired replicates (i.e. non-parametric version of 1-way ANOVA).

Non-parametric tests often lack power and I suspect for this sample size that will be the case. Personally I would try and ignore all this advice :-) - okay I would definately check for normality and investigate avenues towards normality - then I would perform a general linear model (GLM). The reason is that it is a much more powerful method for identifying whether the two gropus are different. Performing a single GLM with two paired groups is not trivial, so I would have to give serious consideration to 1-way ANOVA, or KW-test, which would be a cleaner solution - but could miss a critical part of the data.

"i.e. liver and spleen." after discussing i found biologically it would be more relevant to compare Liver WT and Liver KO and same goes with spleen .Meanwhile i will check the normality of the data which i never used to do because lack of basic stat concept . So to make it more clear since i have sample size uneven in number .so if it passes the normality criteria i can go with unpaired t test ?

What you are saying is you want to do 2 unpaired T-tests, okay fine. I would much prefer to do one test with 2 groups (liver and spleen) using 1-way ANOVA. It might also give you more statistical power. You can always try it with 2 t-tests and 1 one-way ANOVA test. Obviously if you can combine liver and spleen, its a single unpaired T-test.

"would much prefer to do one test with 2 groups (liver and spleen) using 1-way ANOVA. It might also give you more statistical power. You can always try it with 2 t-tests and 1 one-way ANOVA test." thank you for your valuable insight

So far what i know is i have to perform a parametric test .

How did you determine this?

But im not sure if I have to go for a
Paired t-test or Unpaired t-test

If you can't tell the difference between when to do a paired t-test and an unpaired one, you need to stop and talk to someone who can teach you basic statistics.  Throwing statistical tests at numbers when you don't understand them is a recipe for disaster.

well "difference between when to do a paired t-test and an unpaired one" this part i know when to do. But my doubt as i have difference in sample number for WT vs KO cases .Here in this case it would be unpaired which i tried and there is no difference but when i ran it with paired it shows significance. Thats why i would like to know if im doing it right or not

You can't run what you know is the wrong test just because you like the p-value!

"You can't run what you know is the wrong test just because you like the p-value" .perhaps yes


text:
How to convert a .fastq chip-seq file into a wig/bigwig track and how to extract chip-seq peaks?
I have a reference .fasta file and a raw .fastq file with chip-seq data. I am trying to create a bigwig track from the and .fastq and .fasta ref file of the raw signal. Then I would like to do some peak analysis for the chip-seq track.
Can someone help me with googling the write things, I feel really lost while trying to google it and not sure where to start from. I was thinking to do .fastq -> .bam -> wig -> bigwig and then somehow to do peak detection analysis.

Maybe follow some tutorials first? On Biostars many links to helpful tutorials.

Before you proceed, it is usual to do quality control on the .fastq data.  The FASTQC program is a good place to start.
Once you have checked the quality, trimming is the typical next step to remove low-quality data.  There are options such as Trimomatic, trimGalore for this.
Once you are sure you have clean data, your first step is alignment - align the sequences in the .fastq file to the reference sequence in the .fasta file.  bowtie2 is one option for this process.  This gets you to the .bam stage.
These steps will take a little time, some research, and may be somewhat frustrating, but they are the foundation for what you want to do. One option for somewhat automating these steps is to use the Galaxy platform (usegalaxy.org), otherwise these are mostly command-line tools.
The bam to wig and wig to bigwig steps are relatively simple conversions, and Biostars has answers for these.

What about peak calling once I have a bigwig file?

I have heard of the R package BayesPeak, but ChIP-Seq is not really my thing.

macs2 should do the work.


text:
Is there a convention or standard of the highest bits in a sequence logo?
In bioinformatics, a sequence logo is a graphical representation of the sequence conservation of nucleotides (in a strand of DNA/RNA) or amino acids (in protein sequences).
wiki says

The total height of the letters depicts the information content of the
  position, in bits.

and then gives this example

A sequence logo showing the most conserved bases around the initiation
  codon from all human mRNAs (Kozak consensus sequence). Note that the
  initiation codon is not drawn to scale, otherwise the letters AUG
  would each have a height of 2 bits.

it seems that there is a 0.1 bit letter.
is this in the context of 2bits? is 2bits the highest bits in this sequence logo? Is there a convention or standard of the highest bits in a sequence logo?

Sequence logos are usually just visual aids, why do you need the bits? I mean, if you want to get numbers representing the frequency of a residue at a given position, then logos are probably not the right approach.

This is explained in the original paper [1] describing the logos (emphasis mine):

The importance of a particular position in a binding site is more clearly
  and consistently given by the information required to describe the pattern
  there [7, 8]. To choose one symbol or state from two equally likely possibilities
  requires one "bit" of information. For example, to communicate the result of
  a coin-ip to someone requires 1 bit of information because only one yes-no
  question needs to be answered: "Is it heads?". If a position in a binding
  site always contains one base (e.g. G), then we need exactly two bits of
  information since two yes-no questions need to be answered: "Is it A or G?"
  (i.e. is it a purine?) and "Is it A or C?". (If the answers to both questions
  are "no", it must be a T.) Furthermore, if a position contains two bases (e.g.
  sometimes A and sometimes G), only one question suffices since a two out of
  four choice is equivalent to a one out of two choice. Therefore, only one bit
  is needed to describe a position in a binding site that contains only purines,
  but two bits are needed to describe a position that always contains adenine.

So yes, 2bits is the highest value you will see on the y-axis of a DNA sequence logo (proteins are a different story since there are more possibilities). In the one you show, the scale is from 0 to 0.2, but it's the same principle.
However, you should think of logos as a useful visual aid, and not a way of getting rigorous mathematical information about a sequence. So don't put too much weight on the bit values. If you need that sort of data, calculate it from the alignment directly.

[1] Schneider TD, Stephens RM., Sequence logos: a new way to display consensus sequences., Nucleic Acids Res. 1990 Oct 25;18(20):6097-100.


text:
How to reproduce this specific sequence logo from 10 000 human mRNAs?
this is A sequence logo showing the most conserved bases around the initiation codon from 10000 human mRNAs.

I would like to reproduce this sequence logo. How can I get the sequences? which database should I use?
this search is apparently wrong, since I don't know the right way.
https://www.ncbi.nlm.nih.gov/nuccore/?term=10000+human+mRNAs

The figure has no reference to any study, only a website. However on that website no mention of this study. Best thing is to send an email to  weblogo, and ask how they made the figure (I have seen an email address on their site).

If you look at the Appendix from https://www.ncbi.nlm.nih.gov/pmc/articles/PMC306349/, there is a table of sequences that were likely used to construct the sequence logo, at least partially, from the 5' end to the stereotyped ATG subsequence.
It's not in a nice format amenable to use directly with weblogo (like FASTA), but with some programming/elbow grease it might be rendered into that form, which you can then use.
Generally, you might investigate the supplementary materials of these and related papers, to see if you might something directly usable.

Thanks for your answer. By "Appendix", do you mean the "Selected References" part, or somewhere else? I didn't find the Appendix section.

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC306349/?page=9

You might go through the examples here, which include some Kozak sequence data: http://teaching.bioinformatics.dtu.dk/teaching/index.php/ExSeqLogos

If you need to generate that exact logo in your question, then you'll need the sequence data from whoever made it.


text:
Access workdir defined on command line from within Snakefile
Snakemake provides access to a workflow object within a Snakefile. This allows one to, for example, have dynamic programmatic access to the directory containing the Snakefile (via the workflow.basedir attribute). Is there a similar way to access the working directory that is specified on the command line using --directory?
Through a bit of exploration I was able to figure out that the workflow object has a workdir() method, but it appears this is used internally by Snakemake for registering working directories, not for access.
UPDATE: It looks like the working directory is stored internally as the _workdir attribute of the workflow object. However, this attribute is unset (default value None) when I'm creating static/global objects in the Snakefile. Is it impossible to access the working directory prior to building the workflow DAG?

You could query Python for the working directory within the Snakefile, since evidently Snakemake changes the actual working directory of the process to the one specified by --directory.    For example, using a stub Snakefile with just these two lines:
import os
print("Current working directory: " + os.getcwd())

And on the command line:
$ snakemake --directory working-directory
Creating specified working directory working-directory.
Current working directory: /home/jesse/dev/snakemake-wd/working-directory
Building DAG of jobs...
Nothing to be done.
Complete log: /home/jesse/dev/snakemake-wd/working-directory/.snakemake/log/2019-08-02T142009.349338.snakemake.log

Thanks! This is the type of thing that seems super obvious after you've pointed it out, but wasn't obvious initially.

No problem.  I wasn't sure myself if that was how they implemented --directory until I tried it out.  I guess it's an implementation detail but seems unlikely to me it'll change.


text:
Clustering using cellranger
Is it possible to run cell ranger clustering using only a gene x cell matrix? Cellranger reanalyze performs clustering but requires a hd5f file. I don't want to run the whole pipeline (including alignment). I was trying to manually generate a hd5f file to input to cellranger reanalyze. Is this possible?

Hi @guest (nice distinctive username), your question is understandable. I'm not sure whether its possible, but seems reasonable. If you are seeking conversion to h5df, is there any chance you can provide a limited example of the raw matrix output/input?

The direct answer is yes. If you have never worked with hdf5 you can start  here.
But, cellranger outputs several hdf5 files. You can extract the count matrix from the  "/outs/filtered_gene_bc_matrices_h5.h5" file.
They provide a lot of information on how to use the hdf5 format so their their documentation is a good place to start.
In addition: you can also load the data into Seurat, monocle or scanpy and follow their tutorials to analyze the data.


text:
What is the name of this type of figure?
This figure comes from this wiki page.

I googled "gene location figure" and "gene location plot", none gets similar results. I also tried search by image, none of the results is close.
so, what is the name of this type of figure?

This kind of visualization showing the banding patterns of chromosomes are called "chromosome ideograms".
You can use for example IdeoViz in R to generate them.

Such a figure is called ideogram, or chromosome ideogram.


text:
Huge variability in blastn speed between same short sequence job
I have been running a couple of blast searches and I am very suprised and somewhat frustrated by the huge variability in job completion time. The specifics of my blast search are as follows:

Running blast using webbrowser (chrome) via ncbi
19 nucleotide query
database: WGS of one particular sequencing project/bacterial strain (so database is VERY small, in total approx 4.5 MB divided over 60 contigs)
parameters are mostly untouched and are automatically optimized for short queries by blast

In one instance, I started the job more than an hour ago and it is still not done. Knowing from experience that this previously hasn't taken that much time, I started a new job in a different browser tab with exactly the same sequence, database and parameters, and now it finished within 20 seconds. What on earth is happening here?

Welcome to Bioinformatics.SE! I assume you're running BLAST on the command line, and not using a web browser to connect to NCBI's BLAST website. Are you running this on the computer on your desk, or are you connecting to a remote computer? It's possible the machine you're connecting to is experiencing intermittent disk issues.

You mention "tabs" so what service are you using? Please [edit] your question and tell us exactly what you are running and how. What blast? If you're using a website, which one? If not, what blast executable? On what kind of computer?

Thanks for the tips, I edited my question a bit providing more information - running blast via NCBI (in both cases)

First of all, if you're using an online service like NCBI's BLAST, then the time to completion will depend on the load the service is experiencing. Your job is ut in a queue and will be started when it's its turn. If the service is under heavy load, then its turn might come later than you expect. Then, if the service is under heavy load, even after launching your job, it might take longer to finish because it won't have enough resources to run quickly.
Since you say that the same job finished quickly on a new tab, it sounds like either the first job was last in the queue or there was just a bug that caused it not be launched or something like that. We can't really help with things that resolve themselves like that.


text:
Biopython Failure to break out of loop
After qblast, I wanted to print only details from the first 3 BLAST output. Yet, the code below printed the details of every BLAST output. Why can't I break out of the loop to achieve my desired output?
for alignment in blast_records.alignments:
    countResults = 0
    for hsp in alignment.hsps:
        if countResults == maxResults:
            countResults = 0
        break
        else:
            print(alignment.title)
            print(alignment.length/ len(seq)*100)
            print(hsp.expect)
            print(hsp.identities)
            countResults += 1

Where are you setting maxResults?

In my case:maxResults = 3

Your indentation is wrong. Python won't run this. Can you make sure you've posted exactly the code you're running?

Assuming that your actual code is this:
for alignment in blast_records.alignments:
    countResults = 0
    for hsp in alignment.hsps:
        if countResults == maxResults:
            countResults = 0
            break
        else:
            print(alignment.title)
            print(alignment.length/ len(seq)*100)
            print(hsp.expect)
            print(hsp.identities)
            countResults += 1

Here you have an outer loop, which is iterating over all the possible database matches, and an inner loop, iterating over the alignments within that database entry. You are currently limiting the inner loop, but not the outer loop, which you probably want to do. Try this:
# Take the first 3 database hits
for alignment in blast_records.alignments[:3]:
    for hsp in alignment.hsps:
            print(alignment.title)
            print(alignment.length/ len(seq)*100)
            print(hsp.expect)
            print(hsp.identities)

There was an indentation bug (its a Python thing). Personally I'd use,
for alignment in blast_records.alignments:
    countResults = 0
    maxResults = 3
    for hsp in alignment.hsps:
        if not countResults == maxResults:
            print(alignment.title)
            print(alignment.length/ len(seq)*100)
            print(hsp.expect)
            print(hsp.identities)
            countResults += 1
        else:
            countResults = 0 # I'd dedinitely remove this in your original code
            break

The loop gives a bit more control if there is any subsequent filtering, species, percentage identity etc.. There are other ways to write it.
Here's an alternative
for alignment in blast_records.alignments:
    maxResults = 3
    for countResults, hsp in enumerate(alignment.hsps):
        while countResults <= maxResults:
            print(\
            alignment.title, "\n",\
            alignment.length/ len(seq)*100, "\n",\
            hsp.expect, "\n", hsp.identities, "\n"\
            )

BTW the enumerate function keeps track of the counting.
Alternative print ...
           out = """ alignment.title,
           alignment.length/ len(seq)*100,
           hsp.expect,
           hsp.identities """
           print (out)

Let me know if there are bugs

Somehow I posted this on the community wiki, could someone remove it please?

Hi MIchael. I removed the CW from your post, but if something like this happens again, please flag the post for mod attention and send us a custom message explaining what you need. I only saw this comment by chance and could easily have missed it.

Hi @terdon got it, thanks.


text:
PCA on large sparse matrix of single cell RNA-seq
I received a large sc-RNA-seq data matrix, as of the nature of sc-RNA-seq compared to bulk-RNA-seq the data matrix is very sparse.
Moreover, due to the fact that this is single cell data the cell progress stages become more dominant and might be dominant in finding the PC's.
I want to apply a dimensionality reduction method, such as PCA or tSNE, in order to enter it to some working algorithm to find principle curve for pathway analysis.
Ive encountered Bailey paper: "Principal Component Analysis with Noisy and/or Missing Data", but not sure if it will help me in the case of sc-RNA-seq.
Would appreciate some reference link or other help
Thanks

Welcome to BioinfoSE. This is a fairly standard question, and its fairly easy as everything is in Seurat v3. Beyond that you will need to give information on the data set shape and the field (human genetics?).

@MichaelG. I'm not familiar with Seurat v3, is it a package that handle sparse data? regarding the data, its expression dataset of tumors it is size ~(20000, 500) but with a lot of empty entries; sorry if I'm not that clear I'm just getting into the field.

This tutorial would be a good start: https://satijalab.org/seurat/v3.0/pbmc3k_tutorial.html


text:
How can extract the list of genes name from the raw data in GEO?
The format of the raw data that upload in GEO is different, like TXT, CEL or other form. How can I extract the list of genes name from these formats? And then how can I normalize them with R?
I would be really glad if anyone write step by step a complete set of codes for R with an example.

Answer from @benn, converted from comments:
You can find a step by step description in the limma user's guide. Just download the CEL files, and follow the user guide.
[please edit to improve this answer, if possible]


text:
Number of reactions per metabolic pathway
What is the quickest way to calculate the number of chemical reactions present in each pathway for a given organism and which online database such as Kegg, MetaCyc,Reactome is more helpful/reliable for this type of statistics?

Hello, I think "present" is too ambitious here. I guess "known" is the best you can hope for. But I have not looked at the network data for a while, so I am not sure where is the knowledge now.

Agreed. However I need to run a global statistics on the number of "known" reactions per pathway...clearly I expect some bias.


text:
design formula question
I am not sure I am building the proper design formula for the question I want to answer
I have the following samples with three factors; clone, the structure and the condition.
clone structure diabetic
1 07        2D      Dia
2 21        2D     Ctrl
3 23        2D      Dia
4 32        2D     Ctrl
5 34        2D      Dia
6 43        2D     Ctrl
7 07        3D      Dia
8 21        3D     Ctrl
9 23        3D      Dia
10 32        3D     Ctrl
11 34        3D      Dia
12 43        3D     Ctrl

I want to pull the differentially expressed genes for structure (2D vs 3D) blocking by clone, and the one differentially expressed by condition
 (regardless of the structure variable).
I am not sure how to build the design in one call, therefore I have produced the first object to pull 2D vs 3D blocking by clone using the following.
des_structure <- DESeqDataSetFromHTSeqCount(sampleTable = samples,
                                       directory = "../data/htseq_geneCounts/",
                                       design= ~ clone+structure)

for the condition controlled by structure I have used
des_condidtion <- ddsHTSeq_diabetic <- DESeqDataSetFromHTSeqCount(sampleTable = samples,
                                                 directory = "../data/htseq_geneCounts/",
                                                 design= ~ structure+diabetic)

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

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

thank you for the feedback.


text:
How to transform and sort the matrix to make a heatmap showing signatures?
I have a matrix data with cells as rows and samples as columns. Here I am giving the data with dput
dput(data)

structure(c(0, 0, 0.0095, 0.0091, 0, 0.0195, 0, 0.006, 0.0023,
0.0035, 0, 0.0306, 0, 0.0859, 0, 0, 0.0059, 0.0229, 0, 5e-04,
0, 0.0116, 0.0173, 0.0076, 0.0133, 0.0344, 0, 0.1263, 0, 0, 0.0378,
0, 0, 0.0198, 0, 0.0268, 0, 0, 0.0044, 0.0607, 0, 0.0073, 0,
0.1343, 0.0016, 0.0864, 0, 0, 0, 0.0296, 0.0122, 0.0329, 0.0204,
0.0582, 0, 0, 0.003, 0.1011, 0, 0, 0, 0, 0.0124, 0.0013, 0, 0.0025,
0.0076, 0.0498, 0, 0.0103, 2e-04, 0.0562, 0.0673, 0.069, 0, 0.0297,
0, 0.0207, 0.0329, 0.0072, 0.0028, 0.0589, 0.0533, 0.0139, 0,
0.0791, 0.0361, 0.026, 0, 0.0036, 0, 0.0319, 0.0155, 0.0401,
0.064, 0.0743, 0.0026, 0.012, 0, 0.0089, 0.0468, 0.0994, 0, 0.0646,
0, 0.0391, 0.0209, 0.0051, 0.007, 0.0689, 0, 0.0218, 0, 0, 0.0157,
0, 0, 0.0119, 0, 0.0136, 0, 0.0036, 0.0016, 0.0583, 0, 0, 0.0031,
0.2211, 0.6257, 0, 4e-04, 0.173, 0.0557, 0.1104, 0.0288, 0.0474,
0, 0.0078, 0.0914, 0.0289, 0.0293, 0.1067, 0, 0.0121, 0, 0, 0.018,
0, 0.021, 0.039, 0.0146, 0.043, 0.0367, 0.0245, 0, 0, 0, 0, 0,
0.0097, 0, 0.003, 0.0072, 0.017, 0.0061, 0.0453, 0, 0.0587, 0,
0, 0.1009, 0.0231, 0, 0.0279, 0, 0.048, 0.01, 0, 0.0083, 0.0355,
0, 0.1058, 0, 0, 0.0634, 0.0138, 0, 0.0116, 0, 0.0683, 0, 0,
0.0229, 0.064, 0, 0, 0, 0, 0.0222, 0.005, 0, 0.0176, 0, 0.0291,
0.0017, 0.008, 0.003, 0.0693, 0, 0.0156, 0.0038, 0.0577, 0, 0,
0, 0, 0.0521, 0, 0, 0.0017, 0, 0.0218, 0.1336, 0, 0, 0, 0, 0,
0, 0.0398, 0, 0.0099, 9e-04, 0.0046, 0.0055, 0.0361, 0.0382,
0.0168), .Dim = c(14L, 17L), .Dimnames = list(c("Cell1", "Cell2",
"Cell3", "Cell4", "Cell5", "Cell6", "Cell7", "Cell8", "Cell9",
"Cell10", "Cell11", "Cell12", "Cell13", "Cell14"), c("Sample1",
"Sample2", "Sample3", "Sample4", "Sample5", "Sample6", "Sample7",
"Sample8", "Sample9", "Sample10", "Sample11", "Sample12", "Sample13",
"Sample14", "Sample15", "Sample16", "Sample17")))

Annotation of samples is in a dataframe df look likes this:
dput(df)

structure(list(Samples = c("Sample2", "Sample6", "Sample7", "Sample8",
"Sample9", "Sample10", "Sample11", "Sample16", "Sample17", "Sample1",
"Sample3", "Sample4", "Sample5", "Sample12", "Sample13", "Sample14",
"Sample15"), Group = c("A", "A", "A", "A", "A", "A", "A", "A",
"A", "B", "B", "B", "B", "B", "B", "B", "B")), row.names = c(NA,
-17L), class = c("tbl_df", "tbl", "data.frame"))

Using the above data and annotation with below code I made a heatmap:
#Associated libraries
library(ComplexHeatmap)
library(RcolorBrewer)
library(circlize)

#Set annotation
ann <- data.frame(df$Group)
colnames(ann) <- c("Group")
colours <- list("Group"=c("A"="orange","B"="darkgreen"))
colAnn <- HeatmapAnnotation(df=ann, which="col", col=colours, annotation_width=unit(c(1, 4), "cm"), gap=unit(1, "mm"))

myCol <- colorRampPalette(c("#ffffff","#0030dd", "#3259e3", "#7f97ee"))(100)
myBreaks <- seq(0,1, length=100)
hmap <- Heatmap(data, name = "xCell scores", col = colorRamp2(myBreaks, myCol),
                width = unit(100, "mm"),show_row_names = TRUE, show_column_names = FALSE, cluster_rows = FALSE,
                cluster_columns = FALSE, show_column_dend = FALSE, show_row_dend = FALSE,
                row_dend_reorder = FALSE, column_dend_reorder = FALSE, clustering_method_rows = "ward.D2",
                top_annotation=colAnn)
draw(hmap,heatmap_legend_side="right", annotation_legend_side="right")

The heatmap looks like this:

But I want to sort the rows (cell1 to cell14) keeping the cell4 (interested cell) on the top and also keep the samples with values (coloured in blue in heatmap) together and samples with 0 values (white coloured) following that.
It should look something like the below heatmap which is from the Figure 4a in this research paper
As you see in this heatmap all the coloured samples are together.
This must be done with some sorting and using transform, but not able to make it.
Any help is appreciated. thank you in advance.

complex heatmap oncoprint might help but will run it on your data and get back to  you.

The pattern in Figure 4A is created by sorting row 16 (unknown) in each of the 3 groups (Deleterious, Missense, and Wild-type). This new order can then be given to the Heatmap function with argument column_order = vector, where vector is this new order.

Yes, I know this but there are more than 60 cell types. can't make a vector of all types...

Okay, I understand what you mean now, you might want to use a nested sort of let's say the first 10 cell types (depending on how many columns you have in the real data). Nesting can be done by adding the following cell types in the order function.

You are almost there except for a few tricks:
If you want columns to be ordered using the values in Cell 4, you will need to provide a custom column order either with the column_order argument of the Heatmap() function or by providing ordered data. Below is the code to perform the latter. Most of the code is yours, I just tweaked it a little.
library(tibble)
library(dplyr)
library(ComplexHeatmap)
library(RColorBrewer)
library(circlize)

# Since you want to use cell4 as reference
cell4 <- as.data.frame(t(data["Cell4",]))

df <- cbind(cell4,
            Group = as.character(ann$Group))

# ordering by Group and by expression value, minus denotes decreasing order
df <- df[with(df, order(Group,-Cell4)),]

# to make it compatible to your code in your post
df <- rownames_to_column(df, "Samples")

#Set annotation
ann <- data.frame(df$Group)
colnames(ann) <- c("Group")
colours <- list("Group" = c("A" = "orange", "B" = "darkgreen"))
colAnn <-
  HeatmapAnnotation(
    df = ann,
    which = "col",
    col = colours,
    annotation_width = unit(c(1, 4), "cm"),
    gap = unit(1, "mm")
  )

# rownames of the annotation df and the colnames of the data should be in
# the same order
data <- data[, df$Samples]

myCol <-
  colorRampPalette(c("#ffffff", "#0030dd", "#3259e3", "#7f97ee"))(100)
myBreaks <- seq(0, 1, length = 100)
hmap <-
  Heatmap(
    as.matrix(data),
    name = "xCell scores",
    col = colorRamp2(myBreaks, myCol),
    width = unit(100, "mm"),
    show_row_names = TRUE,
    show_column_names = FALSE,
    cluster_rows = TRUE, # you asked your cells to be ordered
    cluster_columns = FALSE, # column order is the same as the provided data
    show_column_dend = FALSE,
    show_row_dend = FALSE,
    row_dend_reorder = FALSE,
    column_dend_reorder = FALSE,
    clustering_method_rows = "ward.D2",
    top_annotation = colAnn,
    column_order = NULL
  )
draw(hmap,
     heatmap_legend_side = "right",
     annotation_legend_side = "right"

)
The resulting heatmap, still not quite what you want:
You will realize that the values for Cell 4 are ordered from high to low for each of the two groups. To make Cell 4 be drawn at the top, you will need to reorder the rows. To do so, you will need to extract the row orders with the row_order() function:
> row_order(hmap)
 [1]  3  2 14  4 12 11  9 10  7  1  5 13  6  8

Now you can re-draw the heatmap after altering the row order somehow. For this sample case, doing it manually.
# promoting Cell 4 to the first position
modified_row_order <- c(4, 3, 2, 14, 12, 11,  9, 10,  7,  1,  5, 13,  6,  8)

hmap <-
  Heatmap(
    as.matrix(data),
    name = "xCell scores",
    col = colorRamp2(myBreaks, myCol),
    width = unit(100, "mm"),
    show_row_names = TRUE,
    show_column_names = FALSE,
    cluster_rows = FALSE, # set to FALSE as a custom row order is supplied
    cluster_columns = FALSE,
    show_column_dend = FALSE,
    show_row_dend = FALSE,
    row_dend_reorder = FALSE,
    column_dend_reorder = FALSE,
    clustering_method_rows = "ward.D2",
    top_annotation = colAnn,
    column_order = NULL,
    row_order = modified_row_order # argument to supply custom row order
  )
draw(hmap,
     heatmap_legend_side = "right",
     annotation_legend_side = "right")

If I have misunderstood your question and you do not want columns to be ordered according to expression values in Cell4, forget about the first part with the ordered annotation data frame.

You may also use the excellent pheatmap package. It is highly customizable, so do not hesitate to look at its help page and examples found therein. The code may look like this:
library(pheatmap)
data<-as.data.frame(data)
df<-data.frame(df, row.names=1)
pheatmap(data, annotation_col = df, cluster_cols = F, color=colorRampPalette(c("white", "steelblue"))(100))


text:
DNA bioinformatic search and methylation analysis by qMSP
Im quoting the methods which was used

We performed a bioinformatics analysis to identify potential CpG
  islands in the genomic sequences of RANKL and OPG genes. The USCS
  Human Genome Browser public database (genome.ucsc.edu), the Primer
  Express software (available from Applied Biosystems), and CpG Island
  Explorer software (cpgie.sourceforge.net) were used to explore the
  nucleo- tide sequences surrounding de transcription start site (TSS)
  of RANKL isoforms I and II (NG_008990.1 Ref Seq Gene), and OPG
  (NG_012202.1 Ref Seq Gene). A (C + G)/total bases ratio

0.5 and CpG observed/CpG expected ratio >0.6 were used as criteria.

Im looking for the region which they have used
This is the paper
"Role of DNA methylation in the regulation of the
RANKL-OPG system in human bone"

I tried using ucsc to check the region but didn't get much. The idea is to find the region to design the primer for the same
How can i get the region from ucsc ,the CpG island around the TSS of the above mentioned genes which are RANKL and OPG
Any suggestion or help would be highly appreciated

It's likely that they accomplished this by using the UCSC Table Browser to:

Find the loci of the two isoforms you mention
Finding CpG islands (Group = Regulation, Track = CpG Islands) near those transcription start sites
Design primers using Applied Biosystem's software

Alternatively, they list the exact primers used in Table 1.
They don't detail exactly how they did it, so if you have specific questions it might be worth it to email the corresponding author

thank you for the suggestion i would try it out and yes it would be better if i drop them a mail


text:
Error using bseqsc
I will be very grateful for any hint on how to overcome the error. I wish to deconvolve my bulk RNA seq data obtained from the lungs of mice using single cell RNA seq data. For practice, I am following this tutorial using my own bulk seq data and the single cell rna seq data from the pancreas used in the tutorial. All seems to go well until when I try to obtain the cell type-specific differential expression by running the following code. From the beginning I thought the presence of factors(e.g. for the treatment conditions or phenoData slot) in my data set was the cause but after getting rid of the factors I still get the same error:
csfit <- bseqsc_csdiff(the_best_ilc_bulk_expr_set[genes, ] ~ phenoData| alpha + beta + ductal + acinar,
                       verbose = 2, nperms = 5000, .rng = 12345)

Then I get this error:
Groups: bleo=5L | bleo_ko_i=5L | bleo_wt_i=5L | Nacl=5L | NA0L
Cell type(s): 'alpha', 'beta', ..., 'acinar' (4 total)
Fitting mode: auto
Data (filtered): 1202 features x 20 samples
Model has factor effect with more than 2 levels: fitting lm interaction model
Fitting model with nonnegative effects
Model with more than 2 groups: switching to version 2
 Fitting linear interaction model ... Error in dimnames(covmat.unscaled) <- list(xnames, xnames) :
  length of 'dimnames' [1] not equal to array extent
In addition: Warning messages:
1: In lsfit(D, G, intercept = FALSE) : 'X' matrix was collinear
2: In sqrt(((n - p) * stddevmat^2 - resids^2/(1 - hatdiag[good]))/(n -  :
  NaNs produced

Here are the details of my code:
# inspired from https://github.com/cozygene/bisque/issues/4

# creating expression set from my bulk rna seq data

library(openxlsx)
library(Biobase)
library(bseqsc)
library(tidyverse)
bseqsc_config('CIBERSORT.R')
data(pancreasMarkers)

# read in the raw data

bulk <- read.xlsx("all_samples.xlsx",
                  sheet = "rawCounts_unrepeated_genes")
head(bulk)

bulk$Gene <- NULL
rownames(bulk) <- toupper(bulk$external_gene_name)
bulk$external_gene_name <- NULL
bulk <- as.matrix(bulk)
head(bulk)

# create expression matrix

featureData <- as.character(rownames(bulk))
featureData <- as(data.frame(featureData,
                             stringsAsFactors = FALSE), "AnnotatedDataFrame")
rownames(featureData) <- rownames(bulk)
phenoData <- c(rep("bleo", 5), rep("Nacl", 5), rep("bleo_wt_i", 5), rep("bleo_ko_i", 5))
phenoData <- as(as.data.frame(phenoData), "AnnotatedDataFrame")
rownames(phenoData) <- colnames(bulk)
the_best_ilc_bulk_expr_set <- ExpressionSet(assayData = bulk,
                                            phenoData = phenoData, featureData = featureData)

# read in the single cell rna seq data and fit the model

eislet <- readRDS('islet-eset.rds')
B <- bseqsc_basis(eislet, pancreasMarkers,
                  clusters = 'cellType', samples = 'sampleID', ct.scale = TRUE)
fit <- bseqsc_proportions(the_best_ilc_bulk_expr_set, B, verbose = TRUE)

pData(the_best_ilc_bulk_expr_set) <- cbind(pData(the_best_ilc_bulk_expr_set), t(coef(fit)))

fit_edger<- fitEdgeR(the_best_ilc_bulk_expr_set, ~phenoData,
                     coef = c("phenoDatableo_ko_i", "phenoDatableo_wt_i", "phenoDataNacl"))

# extended

fit_edger_ext <- fitEdgeR(the_best_ilc_bulk_expr_set, ~ phenoData + beta + ductal + acinar +gamma,
                          coef = c("phenoDatableo_ko_i", "phenoDatableo_wt_i",
                                   "phenoDataNacl", "beta", "ductal",
                                   "acinar", "gamma"))

fit_edger_ext$Symbol <- rownames(fit_edger_ext)

# gather P-values from both models

df_fit_edger <- as.data.frame(fit_edger, stringsAsFactors = FALSE)
df_fit_edger$Symbol <- rownames(df_fit_edger)

req_df_fit_edger <- df_fit_edger["PValue"]
colnames(req_df_fit_edger) <- "Base"
head(req_df_fit_edger)
req_df_fit_edger$Symbol <- rownames(req_df_fit_edger)
df_fit_edger_ext <- as.data.frame(fit_edger_ext, stringsAsFactors = FALSE)
df_fit_edger_ext$Symbol <- rownames(df_fit_edger_ext)

head(df_fit_edger_ext)

req_df_fit_edger_ext <- df_fit_edger_ext["PValue"]
colnames(req_df_fit_edger_ext) <- "Adjusted"
req_df_fit_edger_ext$Symbol <- rownames(req_df_fit_edger_ext)

edger_pvals <- req_df_fit_edger_ext %>%
  inner_join(req_df_fit_edger)
head(edger_pvals)#

rownames(edger_pvals) <- edger_pvals$Symbol

edger_pvals <- mutate(edger_pvals, Regulated = Adjusted <= 0.01 & Adjusted <= Base)

# plot Base vs Adjusted
pvalueScatter(Base ~ Adjusted, edger_pvals, pval.th = 0.01, label.th = 3.5)

# ER genes
genes_ER <- c('HSPA5', 'MAFA', 'HERPUD1', 'DDIT3', 'UCN3', 'NEUROD1')
# Fit on ER stress genes and genes regulated beyond cell type proportion differences
genes <- union(genes_ER, subset(edger_pvals, Regulated)$Symbol)

csfit <- bseqsc_csdiff(the_best_ilc_bulk_expr_set[genes, ] ~ phenoData| alpha + beta + ductal + acinar,
                       verbose = 2, nperms = 5000, .rng = 12345)

Then I have the error:
Groups: bleo=5L | bleo_ko_i=5L | bleo_wt_i=5L | Nacl=5L | NA0L
Cell type(s): 'alpha', 'beta', ..., 'acinar' (4 total)
Fitting mode: auto
Data (filtered): 1202 features x 20 samples
Model has factor effect with more than 2 levels: fitting lm interaction model
Fitting model with nonnegative effects
Model with more than 2 groups: switching to version 2
 Fitting linear interaction model ... Error in dimnames(covmat.unscaled) <- list(xnames, xnames) :
  length of 'dimnames' [1] not equal to array extent
In addition: Warning messages:
1: In lsfit(D, G, intercept = FALSE) : 'X' matrix was collinear
2: In sqrt(((n - p) * stddevmat^2 - resids^2/(1 - hatdiag[good]))/(n -  :
  NaNs produced

Session info:
> sessionInfo()
R version 3.6.0 (2019-04-26)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1

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

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

other attached packages:
 [1] edgeR_3.26.6          limma_3.40.6          preprocessCore_1.46.0 e1071_1.7-2           forcats_0.4.0         stringr_1.4.0
 [7] dplyr_0.8.3           purrr_0.3.2           readr_1.3.1           tidyr_0.8.3           tibble_2.1.3          ggplot2_3.2.0
[13] tidyverse_1.2.1       bseqsc_1.0            csSAM_1.4             Rcpp_1.0.2            Biobase_2.44.0        BiocGenerics_0.30.0
[19] openxlsx_4.1.0.1

loaded via a namespace (and not attached):
 [1] nlme_3.1-140         lubridate_1.7.4      bit64_0.9-7          doParallel_1.0.14    RColorBrewer_1.1-2   httr_1.4.0
 [7] tools_3.6.0          backports_1.1.4      R6_2.4.0             DBI_1.0.0            lazyeval_0.2.2       colorspace_1.4-1
[13] withr_2.1.2          tidyselect_0.2.5     gridExtra_2.3        xbioc_0.1.17         bit_1.1-14           compiler_3.6.0
[19] cli_1.1.0            rvest_0.3.4          xml2_1.2.1           pkgmaker_0.28        scales_1.0.0         NMF_0.22
[25] digest_0.6.20        pkgconfig_2.0.2      bibtex_0.4.2         rlang_0.4.0          readxl_1.3.1         rstudioapi_0.10
[31] RSQLite_2.1.2        generics_0.0.2       jsonlite_1.6         dendextend_1.12.0    zip_2.0.3            magrittr_1.5
[37] Formula_1.2-3        munsell_0.5.0        S4Vectors_0.22.0     viridis_0.5.1        stringi_1.4.3        yaml_2.2.0
[43] plyr_1.8.4           grid_3.6.0           blob_1.2.0           crayon_1.3.4         lattice_0.20-38      haven_2.1.1
[49] splines_3.6.0        hms_0.5.0            locfit_1.5-9.1       zeallot_0.1.0        pillar_1.4.2         rngtools_1.4
[55] reshape2_1.4.3       codetools_0.2-16     stats4_3.6.0         glue_1.3.1           BiocManager_1.30.4   modelr_0.1.4
[61] vctrs_0.2.0          foreach_1.4.7        cellranger_1.1.0     gtable_0.3.0         assertthat_0.2.1     gridBase_0.4-7
[67] xtable_1.8-4         broom_0.5.2          class_7.3-15         viridisLite_0.3.0    iterators_1.0.12     AnnotationDbi_1.46.0
[73] registry_0.5-1       memoise_1.1.0        IRanges_2.18.1       cluster_2.1.0

I have the same problem. Have you solved it? Can you give me some advices?


text:
significant SNPs to annotated candidate genes
I have performed GWAS and got some significant SNPs. How can I get candidates genes in the flanking regions of those SNPs and also their functional annotations (eg. KEGG).

Hi and welcome to BioinfoSE! The amount of information in the post is far too low to answer it. What API you are using? What packages are you familiar with? Any information on biological system is also very helpful. Could also clarify the question, it infers that you already have the genes for the SNPs and annotation of the gene?

hi. I used TASSEL and the R package GAPIT to perform GWAS and got a list of significant SNPs  with info such as chromosome, location (bp), p value. From these SNPs now I need to identify candidate genes from flanking regions of these SNPs and also need to annotate these genes according to their function which will make it easy to identify the potential candidates. So far information from @plat has been helpful and any further help is welcome.

You can use BEDOPS' closest-features tool. You only need to have the SNPs coordinates in .bed format and a .bed file with gene locations. BEDOPS will look for the closest genes upstream and downstream of your SNPs. The command will look something like this:
closest-features yourSNPs.bed yourAnnotation.bed

Then, you can obtain your functional annotations through different tools such as biomart.

OK thanks very much for the help. Will try this one out.


text:
Single sample in group: normal pipeline or Kal's Z test
As stated, which one is better for differential expression analysis? When I say normal pipeline I mean limma-voom, edgeR and DESeq2 pipeline for standard analysis. Kal's z test is mentioned in this paper and I know some pre-made tools like the one on QIAgen uses this as well for single vs single comparison
The situation I have is that I have RNA-Seq data for some patients with rare translocations in leukaemia. There are two of these translocation with almost 10 or even more patients so I have no problem with them. And then there are these translocation which mostly I have only 1, some with 2 or 3 patients, for 7 different translocations.
Our hypothesis is that the two larger group can certainly form their own subtypes, and the others may fall into one of the three cases: (1) belongs to one of the two groups despite of the difference in translocation; (2) belongs to a separate group on their own i.e. complete outlier; (3) belongs to a separate group with other members in the "others" i.e. forming a third subtype. And this can be tested with clustering.
For that, I need to do supervised feature selection to filter out useless genes. I can do that with those with multiple samples, but not those isolated cases. So I thought I can get some DGE analysis and use a certain adj.p and fold change as cut off, or adj.p as ranking matrix, and select those genes for clustering.
So my main questions are:
(1) when using limma/edgeR/DESeq2, should the model design include all the translocation labels in one go or should I just add the one contrast I want to make and build a new model everytime?
(2) Is Kal's z test useful for one vs all DGE analysis and is it better?
(3) Are there other outlier analysis working better under this situation?
Many thanks!

You'll have a bit more power if you use all of the samples at once.
Nothing is useful for this, don't waste your time with it.
You're not looking for outliers, you're looking to see how samples cluster. You can just use correlation for that.


text:
How to download all reactions of the TCA cycle from the RHEA database with the SPARQL endpoint?
I am looking for a way to use the SPARQL endpoint of RHEA to download all reactions of the TCA cycle in CSV or JSON format.


text:
How to get a tabular summary from srspair alignment file?
I checked EMBOSS but there is no tabular output for pairwise alignment.
I can write a script but would like to know if there is already a tool to summarize an alignment in srspair format in the form of a table (merger from the EMBOSS collection gives a tabular summary of the differences but needle doesn't).