text:
Experimental Design for Differential expreression analysis
I have a Normal esophageal Fibroblasts (NOFs) cultured in DMEM media; The same NOF also have been cultured with a tumor sample from a patient named 005 on DMEM media; I have also Cancer Associated Fibroblasts (CAFs) cultured on DMEM, on DMEM + tumor from a patient named as 005 and DMEM + tumor from a patient named 036 something like this picture;

I want to know if NOF and CAF are different in terms of gene expression but I am not sure about design for DESeq2 or Edger and ANOVA.
By @Devon's help I have done so
> mycols
   cell_type co_cultured
G2       NOF          NO
G3       NOF         YES
G4       CAF          NO
G5       CAF         YES
G6       CAF         YES

dds=DESeqDataSetFromMatrix(countData = NOFCAF,colData = mycols, design        =~ cell_type*co_cultured)

Now, I am wondering what would be a suitable contrast for results function in DESeq2
> resultsNames(dds)
[1] "Intercept"                   "cell_type_NOF_vs_CAF"        "co_cultured_YES_vs_NO"
[4] "cell_typeNOF.co_culturedYES"

> res=results(dds)
> res
log2 fold change (MLE): cell typeNOF.co culturedYES
Wald test p-value: cell typeNOF.co culturedYES
DataFrame with 2545 rows and 6 columns

Thanks a lot for correcting me

Do you have replicates of any of these?

Following Devon comment, you don't need biological replicates, technical replicates are worthy. It would help to know how many samples do you have available of each type. In addition, it is important to consider potential batch effects (where all samples extracted and sequenced at once?) and con-founders (are some sample from older patients? Is sex equally distributed)

No, you really need biological replicates.  The conclusions you can draw from technical replicates only is limited.

Thank you, age and sex of patients is the same. In total I have 5 samples (one for each). The sequesing has been done with HTG EdgeSeq Oncology Biomarker Panel for all. Now, I have a matrix of raw read counts with 5 columns and 2545 genes. I want to find the differentially expressed genes between NOFs and CAFs. Any help please?

You can figure out that error from DESeqDataSetFromMatrix yourself, you shouldn't need any help (hint: design is stated incorrectly).

Sorry, will I need the cell_type , co_culture and interaction of them or I just need the interaction of cell_type and co_culture to find the genes different between NOF and CAF?

You don't require any further help on this, you need to fly free a bit more.

You want name=cell_type_NOF_vs_CAF in results().

@swbarnes2 sometimes it is impossible to get biological replicates (with human samples). Two different biopsies from the bowel wouldn't be medically advised, or it is not done initially and then later on it is impossible to have it without adding more complexity to the analysis IMHO. That said, if you can get more replicates (preferably biological and technical) better.

You're really running out of degrees of freedom and the only actual replicate is a single sample from another patient, so please take the results with a huge grain of salt.
You can use a design of ~cellType * coculture where cellType has levels NOF and CAF and coculture has levels yes and no (ideally you'd block by patient, but that's not an option for you).

Sorry, @Devon, this design would return the interaction of NOF and CAF with co_culture or returns gene different between NOF  and CAF by considering the interaction of cell_type and co_culture? Thank you for assigning your time to me

No, this design fits the interaction as well as the main effects, what coefficient you extract is up to you.


text:
Proof of Breakpoint Reversal Sorting Approximation Algorithm
I'm kind of new to bioinformatics and trying to self-study. I'm reading a bioinformatics book: An Introduction to Bioinformatics Algorithms and ran into some problems about understanding the proof of the Breakpoint Reversal Sort algorithm in chapter 5.4.
Specifically, the book states a theorem:

If a permutation $\pi$ contains a decreasing strip, then there is a reversal $\rho$ that decreases the number of breakpoints in $\pi$, that is, breakpoint ($\rho(\pi)$) < breakpoint($\pi$).

This theorem guarantees that the algorithm will terminate and lead to no breakpoint. The book later explains the proof for the theorem:

Let $k$ be the smallest element of decreasing strip, so element $k-1$ would be the ending of an increasing strip. Therefore element $k$ and $k-1$ corresponds to a breakpoint. Then reversing the segment between $k$ and $k-1$ will decrease the number of breakpoints.

I'm quite confused about this proof and don't really understand what the authors mean by reversing the segment between $k$ and $k-1$. Can someone please explain the proof for me?
Here's the lecture slide about reversal sorting that based on the book http://csbio.unc.edu/mcmillan/Media/Lecture07Spring2015.pdf.

Check out this lecture for an easier understanding.

Denote the strips in your permutation as $\langle S_1, S_2, \cdots, S_n \rangle$. Suppose $k-1$ is in segment $S_i$, and $k$ is in segment $S_j$, for some $i, j \in [n], i \neq j$.
Recall by definition of $k$ that $k-1$ must be the maximum and therefore rightmost element in an increasing strip, and $k$ the minimum and therefore also rightmost in a decreasing strip. Consider the two cases:

$i < j$. The increasing strip comes before the decreasing. $k-1$ is on the 'right' end of $S_i$, but $k$ is on the 'right' end of $S_j$ also. In order for $k-1$ and $k$ to be adjacent, reverse the segments $S_{i+1}, \cdots, S_j$.
$j < i$. The decreasing strip comes before the increasing. $k$ is on the 'right' end of $S_j$, but $k-1$ is on the 'right' end of $S_i$ also. Again, we reverse $S_{i+1}, \cdots, S_j$.

So, "reversing the segment between $k$ and $k-1$" means to reverse all the strips in-between the two strips containing $k$ and $k-1$, excluding the 'leftmost' strip, but including the 'rightmost' strip.

source: Stack Exchange

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

sort: votes


text:
Can blat use more than one core/CPU to speed up the alignment?
I am using BLAT to align two versions of the genome of C. elegans. I can see in the Activity Monitor of my Mac Book Pro High Sierra that blat is using 100% of a CPU. However, is this programme able to use more than one core / CPU to speed up the alignment?
I cannot find such option in the help page of the tool here.
I installed blat using conda install -c bioconda blat.
I am using the latest version of blat, i.e. 36.

BLAT can only use one CPU. It is actually not the right tool for full-genome alignment. For "two versions" of the same species, MUMmer and minimap2 are orders of magnitude faster and probably give better alignment.
EDIT (moving comment to answer): OP comments that the purpose of this alignment is for lifting over annotations using the UCSC PSL-based pipeline. So far there are no converters for minimap2's SAM/PAF or MUMmer's delta formats. With the UCSC pipeline, you have to write a converter by yourself. An alternative approach is to use minimap2's liftover. The minimap2 cookbook gives an example:
minimap2 -cx asm5 --cs ecoli_ref.fa ecoli_canu.fa > ecoli_canu.paf
echo -e 'tig00000001\t200000\t300000' | paftools.js liftover ecoli_canu.paf -

Note that this pipeline is not as mature as UCSC liftOver and is not as scalable.

I need an alignment in PSL format for downstream tools needed for lifting over an annotation. Do you know if I can get an alignment in psl format with nucmer or minimap2?

@charlesdarwin LASTZ tool might also be worth a look for this, although it's quite old now...

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


text:
How can I install bioconductor-gviz and use it in jupyter notebook?
I tried to install gviz in a conda environment, but that library seems to be incompatible to python and r. I tried to setup a clean environment using
conda create -n r -c conda-forge r-essentials jupyter

and then add the library with:
source activate r
conda install -c bioconda bioconductor-gviz

getting

UnsatisfiableError: The following specifications were found to be in conflict:
    - atk

I remove atk, but now I get:

UnsatisfiableError: The following specifications were found to be in conflict:
    - bioconductor-gviz
    - r-bindr

Does anyone manage to use gviz from within an jupyter notebook?
I also tried to install gviz from a running R-notebook with:
install.package('gviz')

Warning in install.packages :
    package ‘gviz’ is not available (for R version 3.4.3)

Same when I try 'bioconductor-gviz'.

The package is from Bioconductor, from the version you are using you should use:BiocInstall::biocLite("gviz") to be able to install it from R terminal

Not sure what that means.

That to install a package from Bioconductor you need to use a different command (not install.package), but the one I provided

It does not seem to be an accepted command though.Error in loadNamespace(name) : there is no package called ‘BiocInstall’

Did you try following the installation instructions?
Try this in R:
source("https://bioconductor.org/biocLite.R")
BiocInstaller::biocLite(c("Gviz"))

Nope, I did not know, each package has its own way of installation. But the first line here was what I was missing.

All Bioconductor packages are installed the same way. Though in general, reading the software documentation before use is always a good idea. They've changed the command from my answer a bit for R 3.5 and later (see Alex's answer), but for R < 3.5, this is the general way to install any Biocondoctor package.

# Install the BiocManager package first
install.packages("BiocManager")
library("BiocManager")

# Then install the Gviz package
BiocManager::install("Gviz")

Use BiocManager::install()
# Code
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("Gviz", version = "3.8")

However, If R < 3.5.0 , please use the following for installing Bioconductor packages.
BiocInstaller::biocLite("Gviz")


text:
Circos plot with mutation data
This was my earlier question but now I would like to see patient sample mutation which i have tried and im getting an output but it looks really dumb
So here are my data files and conf file files
In my file which im trying to plot is gene having mutation in terms of patient which ranges from 1  to 53 patient at max .
NPM1 is having the highest mutation with 53 patient in the plot what i tried is to set a threshold of between > 10 and < 3 to set the color label ,which labels the number it contains.

How can i add the gene names to the plot ?
Is there a way where i can show the genes with only 1 mutation are linked together ,similarly for 2 linked together so on and so forth .

I tried with the links but couldn;t make it work.
Any suggestion or help would be highly appreciated .


text:
Using multidimensional scaling to visualize protein sequences by functionality
I have a multiple sequence alignment (MSA) of protein sequences, with which I have performed multidimensional scaling to visualize their clustering.  Using R, I used the seqinr package to create an identity matrix (using dist.alignment).  Is there a similar method that is more "generalizable"?  Specifically, I'm interested in classifying the amino acids by some functional characteristics, such as polarity, instead of having the alignment list the amino acids by their identity (A, G, I, etc.).
Here's a very short set of fake input data, with columns 2-4 representing amino acid residue positions:
a <- c("Strain 1", "Strain 2", "Strain 3", "Strain 4")
b <- c("Polar_hydrox", "Nonpolar", "Nonpolar", "Polar_hydrox")
c <- c("Nonpolar", "Nonpolar", "Nonpolar", "Nonpolar")
d <- c("Nonpolar_cyclic", "Nonpolar_cyclic", "Polar_amide", "Polar_amide")
e <- c("Polar_amide", "Acidic", "Polar_amide", "Acidic")
x <- data.frame(a, b, c, d, e)

Is what I'm asking a feasible/correct use of MDS?  Or should I be looking into a different method entirely?  I attempted to use Gower's formula to create a dissimilarity matrix on the original dataset (i.e. using the actual amino acid identities instead of the functional groupings), but compared to what I got using dist.alignment, the MDS plots look different.
Edit: Gower's formula and the dist.alignment function are obviously calculating the distance matrices differently, but I don't know the details on the different methods.  Would using Gower's formula be a good approach for calculating distance matrices for the example dataset?
Edit: clarified functions used for calculating distance matrix.

If you use different distances metrics it is expected to get different MDS plots. What do you mean by more generalizable? What is your biological goal?

Thanks for your response.  My goal is to see whether sequences appear more related to one another when considering what the functionality of their amino composition is, versus the raw amino acid composition.  I'm interested in the phenomena that not every kind of amino acid change will affect immune escape equally.

The question is then why do you settle for Gower's formula and not other distance metrics like the Manhattan distance? It is as good as any other distance metric as far as I can tell.

Some of the first MSA analysis, eg of codon usage of bacteria were performed using matrix factorization (see old papers from Des Higgins). His group also published discriminative methods for finding functional residues using supervised correspondence analysis

Wallace IM1, Higgins DG BMC Bioinformatics. 2007 Apr 23;8:135.
  Supervised multivariate analysis of sequence groups to identify
  specificity determining residues.

The method is available in the R package ade4 and Bioconductor package made4
Also have a look at the phylogenetic R package adephylo and several related packages

I personally don't think this will work well for some variables, primarily because you are not mapping conformational epitopes, where amino acids will cluster in physical space by their charge. In particular, the charge is critical to retaining this conformational structure.
However, hydrophobic/hydrophilic characters are far more amenable to linear style epitopes that this analysis is primed towards.
I don't recommend a multi-dimensional array for this analysis, because you have loads of information, but rather a 'dictionary' (Python), hash (Perl) (or 'object' [any]) or XML, which retains all the information for a given sequence and then strips this in preparation for your multivariate statistic.


text:
Genomics Statistics Problem
I have used a python script to identify target sequences in a DNA sequence file. There are two classes of sequence: coding and non-coding.
I have identified 728 sequences of interest. 597 of these fall into the coding regions and 131 of these fall into the non-coding regions. This is the equivalent of 18% non-coding. The total non-coding region in the sequence file is 13%.
Is there a statistical tool to demonstrate the python script identified target sequences in a non-random fashion?
If the script identified sequences that were randomly distributed then 13% of them would have been found in the non-coding region, from a total of 728 this seems like it should be reliable.

It sounds like you're fishing for a statistic to prove a point. I wouldn't try to over-analyze your data. You might just consider stating that 18% of the sequences are in non-coding regions which is greater than the expected null hypothesis of 13% of the sequences being in non-coding regions.

One could ask how big are those sequences of interest (SOI) and if there are some of them that span between coding and non-coding regions. Also the distribution of SOI could depend on TADs or other higher structures of your sequence (if it is a crhomosome size sequence). To model the distribution of SOI would require knowing a lot more about the SOI and the sequence you are working with and without a null model it cannot be tested if it is randomly following this model or not

Hi again @conchoecia. Thanks for the feedback. With this point, however, we find ourselves in disagreement. There is a chance that if my results are due to random chance there may well be an uneven distribution. Say if my non-coding hits were 14% it would be obvious this would not be reliable evidence in support of my hypothesis. Furthermore, I have no problem with obtaining a negative result but however I interpret my results it needs to be statistically reliable, there is a chance that 18% is not significantly different to 13% and therefore cannot be used in favour of my hypothesis.

@llrs Hi. The sequences of interest are 29 bp (promotors) and I am searching chromosome FASTA. There are a some overlaps, but I have not counted them as non-coding for simplicity. The non-coding regions were defined by a different study (https://doi.org/10.1186/s13068-017-0742-z) Table 1.

Thanks everyone, question answered on cross-validated.
https://stats.stackexchange.com/questions/386173/genomics-stats-problem/386177#386177

If your null hypothesis is that the event of being selected by your program is statistically independent from the event of being non-coding, you can use Fisher's exact test to reject the null hypothesis.
Suppose you have 10,000 (I'm making that up but presumably you have the real number) total sequences, 8,700 coding and 1,300 noncoding. Given that your program selects 728 sequences, you would expect around 13%*728 ~= 95 noncoding sequences. But your program picked 131 noncoding sequences and you want to know if 131 is far enough away from 95 to conclude that the difference is not just due to chance.
The probability of selecting 131 or more of these noncoding sequences (given that we are choosing 728 total sequences) is given by the so-called hypergeometric distribution. The two-sided test will also incorporate the probably of getting an equally unlikely low number (in this case, somewhere around 57) or less.
In R for example, you can calculate this using the 2-by-2 contingency table as follows:
fisher.test(rbind(c(597, 131), c(8103, 1169)))

which gives:
p-value = 7.399e-05

So in this situation (and assuming the probability model described), there's around a 0.0074% chance of observing a ratio as "extreme" as the one you observed, which indicates it's unlikely that the two events are independent.
Note: As others have mentioned, a binomial test could be appropriate too. If your number of total sequences is very large, then the two results will converge. But if it isn't, then Fisher's exact test will be more appropriate (essentially since you're sampling without replacement).

A Poisson or binomial would suffice.
I'm not at all sure the difference will significantly deviate from randomness.
Perl: 95% confidence interval of a non-parametric bootstrap for e.g. 1000 replicates is an alternative approach. I'd be more comfortable coding that solution in Perl than I would in Python, so is unlikely to be of much use here. Anyway, I would use a 2-dimensional array for your total data-set, the second dimension being coding/non-coding, I would then use the inbuilt random-number generator and multiply by the size of the array to randomly pluck out each element. 95% CI of the mean works by ranking the data and removing the top and bottom 2.5% and taking the interval for a single replicate (which is repeated 999 times more).
There will be a deep learning approach to this data, which would be easier to implement in Python, but I'd need to think about this because it is not simple classification.
Again my warning is that given the sample size, the result at best might be borderline. There are likely to be lots of interesting in your data, but coding/non-coding may not be one of them.

Hi, Michael. Thank you. Yes, I am trying binomial analysis now, don't worry about whether this results in a null hypothesis. I am simply trying to rule out a possible explanation for some results I obtained in the wet-lab. Whatever result I get is fine by me I am just happy to show I have properly investigated possible explanations.

Good to know. I'm definitely advocating DL for molecular solutions, they are very, very cool.

Haha, I'm interested in DL however my data science, coding and statistics skills need a lot of work before I can use such a tool. But I read a lot about the applications and I'm keen to get to the point where I can use such tools.


text:
R package installation in r studio
I have written a package for my personal use and my other colleagues. I have uploaded the files and source package to different websites. I am able to install the source package in my computer. But my colleagues cannot install it either by github.com (https://github.com/furqan915/Epi-gene) or by sourceforge.net (https://sourceforge.net/projects/epi-gene/).
Common error is package cannot be installed and it does not add to user library.
I am unable to find out the problem.
Kindly help me out.

I'm voting to close this question as off-topic because it is about general programming. If you provide the errors you find and the steps you have done to find the cause of them it will be easier for people on Stack Overflow or other places to help you

The package installs without issues on my system. You’ll need to provide a lot more details about this error for people to be able to help you.

Can you use devtools::install_github() giving your github username/reposname
If not, you need to check that your R package is passing R CMD check or devtools::check(). See http://r-pkgs.had.co.nz/check.html and http://kbroman.org/pkg_primer/pages/check.html for info

I have already finished running the RCMD check. I got 2 notes but it finished in good manner.

A quick workaround is to install the package directly from source:
install.packages("~/Download/myPackage.tar.gz", repos = NULL, type="source")

You should be able to pass also .zip and URLs (other than CRAN).


text:
convert supplementary reads to primary in sam or bam
I have a problem with a tool, that possibly ignores supplementary reads. I want to find out a bit how this tool works by converting all supplementary reads to primary reads, and then change the names of those reads (because I guess that primary reads with the same name are ignored?). If the tool produces something, instead of nothing (like now), then I know for sure its the supplementary read flag that could be the problem.
I use SV detection tools. The longer the reads are that I produce (like simulating improvement in the long-read sequencing technology), far less SVs gets detected, while my hypothesis is that longer reads should be able to detect more. Then I found that the alignment files look clean (sam/bam), but with far more supplementary tags. I thought I could test this by changing the flags.
I am searching the internet, but changing read flags is not something they seem to do, only things like omitting reads with certain flags from the file. Maybe I am looking with the wrong search terms, in that case I am sorry.
I guess I need to change column 1 and 2 of the sam file (QNAME and FLAG). But I was also wondering if this is going to work at all. Maybe someone with more knowledge of alignment files could help me with this problem? suggestions how to change the sam/bam file to check if supplementary reads are ignored?
I think this, because I made on purpose a very unrealistic alignment. A 15Mb region against a 15Mb region with differences (x10 times to " simulate" 10 reads / 10x coverage). The alignment file looks quite good actually in IGV, but clicking on any read reveals that there are of course MANY supplementary reads. I also have a "normal" alignment file, where I also want to check if changing the supplementary tag could possibly help me with this tool.


text:
Extracting genes differentially expressed by Wilcox test
I have 2 cell types (NOF and CAF) cultures individually and also have been co-cultured with tumour like this picture.
https://ibb.co/x8Ty1Wz
If this is the log cpm of these 4 samples
> head(onco[,1:4])
             A10       A11        A12        A2
A2M    11.755100 12.145143 12.9787977 12.407221
ABCB11  2.955457  2.278363  0.1229266  1.493507
ABCC2   6.711643  6.974025  0.1229266  7.628936
ABCC6   6.443439  6.094670  6.1389496  7.483394
ABCF1   9.995098 10.045613  9.6627424 10.314842
ABCG2   7.125483  7.448127  6.8669702  8.492895
>

I have used log of count per million (cpm) for Kruskal-Wallis rank sum test that says these samples are significantly difference like below;
> NOF1=c(NOF1$`onco$A10`)
> NOF2=c(NOF2$`onco$A11`)
> CAF1=c(CAF1$`onco$A12`)
> CAF2=c(CAF2$`onco$A2`)
> my_data <- data.frame(
+     group = rep(c("NOF1", "NOF2","CAF1","CAF2"), each = 719),
+     cell_type = c(NOF1,NOF2,  CAF1,CAF2)
+ )

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

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

    Kruskal-Wallis rank sum test

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

>

I then did Pairwise comparisons using Wilcoxon rank sum test to see which samples are different with each other.
> pairwise.wilcox.test(my_data$cell_type, my_data$group,
+                      p.adjust.method = "BH")

    Pairwise comparisons using Wilcoxon rank sum test

data:  my_data$cell_type and my_data$group

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

P value adjustment method: BH
>

My Question is; how I can extract the genes behind the difference between each pair of samples? I mean something like differential expression by non parametric test. also I am not sure if I am right in using log cpm values for these test
If possible give me a hand to correct myself and get intuition
Thanks

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


text:
"Sequence Duplication Levels" module still fails after pre-processing Illumina data
I want to ask about why the sequence duplication levels are high after I trimmed by using Trimmomatic? I am using the following Trimmomatic operations: HEADCROP = 19 TRAILING = 20 MINLEN = 66.
How can i solve this problem? Thank You.

Why do you think this is a problem to begin with?

I though the cross sign (X) means some kind of error and should be eliminated with certain pre-processing technique (like trimming to solve adapter content error)? I am new to Illumina, thank you for advising :)

What type of sequencing libraries are these?  Whole genome, RNA Seq, whole exome, targeted seqeuencing?  What genome? also how complex is the library you sequencing.

To answer your direct question, there are a few reasons why there might be high levels of sequence duplication. From the FastQC help:

The underlying assumption of this module is of a diverse unenriched library. Any deviation from this assumption will naturally generate duplicates and can lead to warnings or errors from this module.

As @DevonRyan mentioned, with certain sequencing protocols such as RNA-Seq, two sequence reads at exactly the same location aren't that uncommon. This isn't a problem with RNA-Seq data, or with Trimmomatic, or with FastQC. It's just that this kind of data violates the assumption, and therefore should be ignored in those circumstances.
PCR duplicates are another possible cause. PCR duplicates can give the false impression of high coverage at a particular locus when in fact it's just a single observed read that has been duplicated many times (see here for more details). PCR duplicates can usually be detected and removed if your analysis involves mapping to a reference genome. But whether this is actually a problem you need to fix depends on what type of data you have and what types of analysis you want to do.
Large numbers of adapter dimers or rRNA may be present in your sample.

But I think it's also important to address how quality control (QC) is run. It can be tempting to run and re-run QC tools like Trimmomatic until all errors go away, but to be blunt these tools cannot think for you. For example, it's possible to get rid of most adapters by aggressively cropping/trimming both ends of each read, but you'll likely throw away a lot of good data that way. You may want to look into Trimmomatic's ILLUMINACLIP operation. It's also may be tempting to crop/trim reads aggressively if there are compositional biases near the beginning or end of the read. In fact, random hexamer priming can cause the Per-base Sequence Content module to fail on almost any RNA-Seq sample. Again, from the FastQC help (emphasis mine):

Libraries produced by priming using random hexamers (including nearly all RNA-Seq libraries) and those which were fragmented using transposases inherit an intrinsic bias in the positions at which reads start. This bias does not concern an absolute sequence, but instead provides enrichement of a number of different K-mers at the 5' end of the reads. Whilst this is a true technical bias, it isn't something which can be corrected by trimming and in most cases doesn't seem to adversely affect the downstream analysis. It will however produce a warning or error in this module.

So in other words, it's best to determine what can cause each FastQC module to fail, investigate whether this is actually a problem for your data set (referring to documentation as needed), and make a deliberate QC plan that addresses the issues that need attention.

FastQC assumes that all samples are for whole genome sequencing and will flag them as failed if they differ too much from that assumption. This will, for example, cause essentially all RNA-seq, ChIP-seq, and ATAC-seq samples to fail in one module or another. This is not any cause for concern and is completely expected. Primarily concern yourself with whether all of your samples are similar in their metrics.


text:
how to write italic script in Rstudio
I have got the problem on how to write the italic script in RStudio. I used the script in the screenshot and got the error "Error in grobs[[i]]: subscript out of bounds."

Please replace the image with the actual text.

Have you ensured that all of the elements of p are grobs (things plottable by ggplot)? The error suggests that that's not that case, so try just iterating over 1:4 instead.

Seriously. It's even easier to paste the text!

I'm voting to close this question as off-topic because formatting text as italics is not a bioinformatics problem

The title is misleading as this error doesn't have anything to do with making the font italic. In Seurat v2, FeaturePlot does not return a ggplot2 object by default, so p in your case is NULL. You need to set do.return to TRUE in the FeaturePlot call.
You should instead do:
library(Seurat)

p <- FeaturePlot(pbmc_small, head(pbmc_small@var.genes), do.return = TRUE)

for(i in 1:length(p)) {
  p[[i]] <- p[[i]] + theme(plot.title = element_text(face = 'italic'))
}

cowplot::plot_grid(plotlist = p)


text:
Highly heterozygous reads mapping
I have short (67bp) Hi-C reads from a highly heterozygous organism (~15% between-haplotype SNP divergence) and I have both reference haplotypes.
I wanted to try and compare different haplotyping softwares for Hi-C reads using these reads as benchmarking dataset. When mapping the reads separately to each haplotype, I get good mapping statistics. When I map the reads to a single haplotype with all the heterozygous SNPs masked (into N's), I get  very poor mapping rates.
I would like to be able to map the reads when the real haplotypes are unknown (reference is a mix of haplotypes).
I use minimap2 to map the reads with the sr preset. I tried decreasing the mismatch penalty (-B) to 1 and increasing the --score-N value, but this had no effect.
As shown in the attached IGV screenshot, coverage drops to 0 when the SNP density increases.
Is it feasible to map reads with such a high heterozygosity on a single (masked) reference ? Should I use another tool ?

Could you describe your reads a bit more? DNA/mRNA? SE/PE? 50bp/100bp/etc? If I were to make a guess, I'd say you're having issues with properly seeding your reads for alignment. Also, it might be helpful to know why you need to map masked reads when mapping separately works fine.

Thanks for the feedback ! I updated the question to add more background and infos.

Just to confirm it is heterozygous and not highly repetitive?

Yup, just highly heterozygous, it is a hybrid yeast

With this levels of heterozygosity you can just assemble everything and get the haplotypes separated. I am not sure if this is an optimal benchmarking dataset.

You're right, I'll be using another dataset for this, but I was still curious about the feasibility :)

I believe you may be able to map your reads, but I don't know how to do that with minimap2.
I recommend running gsnap, which is more SNP tolerant and provides a number of parameters which are likely to help.
For instance, I believe that most aligners will treat 'N' characters as mismatches when you align. GSNAP has parameters to account for this.
--query-unk-mismatch=INT       Whether to count unknown (N) characters in the query as a mismatch
                               (0=no (default), 1=yes)
--genome-unk-mismatch=INT      Whether to count unknown (N) characters in the genome as a mismatch
                               (0=no, 1=yes (default))

It also has a mismatch parameter similar to the one you described for minimap2.
-m, --max-mismatches=FLOAT     Maximum number of mismatches allowed

Try running with the 'genome-unk-mismatch' parameter above (with the masked reference). That may be your best bet. There are also other parameters that may help, but this should be a good start.


text:
making my own population allele frequency table for input to IADMIX from Gnomad data (population admixture)?
Hi I'm a first time user of IADMIX. I tested on one known Finnish sample from the  1000Genome project using the softwares provided frequency table hapmap3.8populations.hg19.freqs and the prediction is correct.
I'm also want to do the same using population frequencies from the gnomAD database but my prediction is wrong - comes out to be African.
My steps are: Create my own population frequency table using the gnomad database from http://gnomad.broadinstitute.org/downloads For each VCF file i do the following
bcftools view -Ou -I -q 0.05  -v snps -k -m2 -M2 -f PASS gnomad.genomes.r2.0.1.sites.22.vcf.gz | bcftools query -f '%CHROM\t%POS\t%ID\t%REF\t%ALT\t%INFO/AF_nfe\t%INFO/AF_afr\t%INFO/AF_amr\t%INFO/AF_eas\t%INFO/AF_asj\t%INFO/AF_fin\t%INFO/AF_oth\n'  -o iadmix-22.freq.bcf

which gives me tabulated data like this containing only "known" variants with RSID != "." , min threshold for total AF >= 0.05 , only bi-allelic SNPs that also have FILTER==PASS

7 16692 rs370739206 G C 1 . . . 1 . 1  7 16719 rs151261124 G A 1 1 1 .
  1 . 1  7 20460 rs1466291194 T C 0.0381589 0.0630435 0.191919 0.0357143
  0.0333333 0.144909 0.121212  7 20471 rs370790239 C T 0.631696 0.716216 0.766667 0.666667 . 0.661765 0.625

Then I have a small R script to - exclude sites that have any population groups AF values == "." - include sites that has mean of all AFs for the snps must be >= 0.05 (to include only common variants)
Then I tried running IADMIX using this frequency table on a know finnish 1KG sample, but the prediction is incorrect
python runancestry.py --freq MY.FREQUENCY.FILE --plink SAMPLE_file --out testoutput --path full/path/

snps in file 455423 pops 7 AF_NFE:0.167 AF_AFR:0.189 AF_AMR:0.133 AF_EAS:0.188 AF_ASJ:0.238 AF_FIN:0.037 AF_OTH:0.047 initial sol 1082621.766160 RUNNING THE L-BFGS-B CODE

    * * *

Machine precision = 2.220446e-16 ... ...

final maxval -1044248.568804 ADMIX_PROP AF_NFE:0.0000 AF_AFR:0.5822 AF_AMR:0.1137 AF_EAS:0.1797 AF_ASJ:0.1244 AF_FIN:0.0000 AF_OTH:0.0000 AF_AFR:0.5822 AF_AMR:0.1137 AF_EAS:0.1797 AF_ASJ:0.1244 FINAL_NZ_PROPS

Any one know whats happening here?
For reference the output using the frequency table provided by IADMIX is as follows

snps in file 106250 pops 8 YRI:0.015 CHB:0.003 CHD:0.275 TSI:0.207
  MKK:0.101 LWK:0.169 CEU:0.194 JPT:0.037  initial sol 76610.617803
  RUNNING THE L-BFGS-B CODE
    * * *

Machine precision = 2.220446e-16 final maxval -72537.379103 ADMIX_PROP
  YRI:0.0000 CHB:0.0000 CHD:0.0000 TSI:0.1031 MKK:0.0000 LWK:0.0000
  CEU:0.8190 JPT:0.0779 TSI:0.1031 CEU:0.8190 JPT:0.0779 FINAL_NZ_PROPS

so the CEU proportion is 82% , acceptable for this finnish sample, right?

I’m voting to close this question because it seems to lack important detail. There's mention of a custom R script (not included) that is used for filtering SNP frequencies, and it's not clear where the rationale behind doing this filtering comes from.

How did you solve this issue? I am having the same exact problem as you


text:
amplicon sequencing set up
I have .fastq files from a targeted amplicon sequencing by Illumina MiSeq run. I want to get the target sequences to perform Blast search on them. What is the pipeline to use for processing amplicon sequencing data? I also have 20% PhiX.

Usually you have the targeted regions based on the library/ panel design prior to sequencing.  Do you not have that information?  It is standard Illumina processing to remove the PhiX reads before creating the FASTQs

This is crossposted on biostars: https://www.biostars.org/p/357824/

Do you want the reads that contain the amplicon? If you know the primer sequences you could do this pretty easily


text:
Drug with as many targets as possible through a system biology graph?
I have 4 pandas data frames containing the following data:

affinity_score1: affinity between each couple (Drug-Target) measured using  tool1
affinity_score2: affinity between each couple (Drug-Target) measured using tool2

The two tools (tool1, tool2) give measure in different scales and are not 100% accurate.

Drug_similarity: a similarity matrix between drugs
Target_similarity: a similarity matrix between targets.

Here also the tool used to measure the similarity is not 100% accurate.
import pandas as pd

affinity_score1 = [{'Drug0': -24.5, 'Drug1': -17.1, 'Drug2': -13.1},
           {'Drug0': -20.0, 'Drug1': -21.6, 'Drug2': -18.1},
           {'Drug0': -19.0, 'Drug1': -15.4, 'Drug2': -20.0}]

affinity_score2 = [{'Drug0': 1.5, 'Drug1': 0.8, 'Drug2': 0.4},
           {'Drug0': 1.2, 'Drug1': 1.3, 'Drug2': 0.9},
           {'Drug0': 0.8, 'Drug1': 0.7, 'Drug2': 0.9}]
affinity_score1_df = pd.DataFrame(affinity_score1).rename(index={0:'TargetX',1:'TargetY', 2:'TargetZ'})
affinity_score2_df = pd.DataFrame(affinity_score2).rename(index={0:'TargetX',1:'TargetY', 2:'TargetZ'})
Drug_similarity = [{'Drug0': 1, 'Drug1': 0.8, 'Drug2': 0.2},
           {'Drug0': 0.8, 'Drug1': 1, 'Drug2': 0.9},
           {'Drug0': 0.2, 'Drug1': 0.9, 'Drug2': 1}]
Target_similarity = [{'TargetX': 1, 'TargetY': 0.8, 'TargetZ': 0.2},
           {'TargetX': 0.8, 'TargetY': 1, 'TargetZ': 0.9},
           {'TargetX': 0.2, 'TargetY': 0.9, 'TargetZ': 1}]
Drug_similarity = pd.DataFrame(Drug_similarity).rename(index={0:'Drug0',1:'Drug1', 2:'Drug2'})
Target_similarity = pd.DataFrame(Target_similarity).rename(index={0:'TargetX',1:'TargetY', 2:'TargetZ'})

I am trying to find the Drug who can have as many Targets as possible.
So far, a simple approach could be to select the Drug having the highest average score in both affinity scoring tools.
This approach does not use the information in the similarity matrices.
Given that if a Drug1 is successful with TargetX, and TargetX is similar to  TargetY then Drug1 can also easily be successful with TargetY. Considering the inaccuracy of the scoring tools, I think including this information could improve the accuracy of the answer, eliminating some false positives.
Also, if a Drug from a cluster of similar drugs (information that can be derived from the drugs similarity matrix) is the only one predicted to have good affinities with some targets, then this is more likely to be a false positive.
But I don't know how to efficiently integrate all the information. Scipy may help with clustering the two similarity matrices.
Another way of analysis is to set affinity and similarity thresholds, in which case each data frame could be transformed into a network:

Each node is a drug or a target.
Two types of edges: similarity (between drugs or targets) and affinity (drug-target)
The connected nodes have passed the thresholds.
The edges have weights: affinity or similarity score.

The problem then becomes finding the most connected node (a drug node) with the highest weight values (scores on tool1 and tool2). Connected nodes through similarity may then be used to validate the number of targets for each drug.
It seems a more complicated approach but more rigorous with most likely higher positive predictive value.
Any python package to implement such an approach?
And also a visualization tool.
The original data has 38 Targets with 806 drugs.

I don't know any python (or R package) to do this), but this is how I would do it based on these matrices:

Do a multidimensional scaling (MDS) plot of the target_similarity matrix.  Explore the visualization:
Find which drug targets the drug closest to the 0,0 point in the 2D space

However this doesn't say anything about the affinity between the drug and the target.
If you wanted to add this then you could weight how much close is the distance between the drug and the target vs how close are they to the center of the MDS (a simple division could work if the distance is scaled between 0 and 1 and averaged).
Note that depending on the original information used to calculate the similarities and the affinities it could be done in other ways

I like the MDS idea even though I don't fully understand the 2nd part.
I can do a MDS on target similarities and another one on drug similarities.
Then create a di-graph like here, part C. Drugs and targets will be connected through edges based on their affinity scores. An affinity threshold will define which drug is connected to which target, and an assigned weight will give the affinity score.

No, I didn't say anything about a MDS on drugs similarities. The center of the MDS of the targets is the closest point to all the targets. That's why you need to find the drug that aims to the target closer to the center of the MDS.


text:
How to plot genomic.fna fasta file in R using Gviz?
I downloaded the genome of Staph aureus as DNA sequence from NCBI. I would like to visualize it using Gviz. Is Gviz the right tool to visualize genomes? What are the necessary steps from the sequence file to an object that can be used in Gviz? I searched the user guide and the reference manual. It is stated that it is possible to use fasta file, but not how to use them.
We will generate lots of genomes that we will have to annotate and compare. S aureus is not S aureus. There is a spectrum of phenotypes, they all have a different genomes. And we will also not be limited to SA. The question here is not how to annotate the genome, but rather the steps necessary to visualize something and get going. It is my first steps with genomic data, I don't even know what makes sense.
Here is the head of the file:
>NC_007795.1 Staphylococcus aureus subsp. aureus NCTC 8325 chromosome, complete genome
CGATTAAAGATAGAAATACACGATGCGAGCAATCAAATTTCATAACATCACCATGAGTTTGGTCCGAAGCATGAGTGTTT
ACAATGTTCGAACACCTTATACAGTTCTTATACATACTTTATAAATTATTTCCCAAACTGTTTTGATACACTCACTAACA
GATACTCTATAGAAGGAAAAGTTATCCACTTATGCACATTTATAGTTTTCAGAATTGTGGATAATTAGAAATTACACACA
AAGTTATACTATTTTTAGCAACATATTCACAGGTATTTGACATATAGAGAACTGAAAAAGTATAATTGTGTGGATAAGTC
GTCCAACTCATGATTTTATAAGGATTTATTTATTGATTTTTACATAAAAATACTGTGCATAACTAATAAGCAAGATAAAG
TTATCCACCGATTGTTATTAACTTGTGGATAATTATTAACATGGTGTGTTTAGAAGTTATCCACGGCTGTTATTTTTGTG
TATAACTTAAAAATTTAAGAAAGATGGAGTAAATTTATGTCGGAAAAAGAAATTTGGGAAAAAGTGCTTGAAATTGCTCA
AGAAAAATTATCAGCTGTAAGTTACTCAACTTTCCTAAAAGATACTGAGCTTTACACGATTAAAGATGGTGAAGCTATCG
TATTATCGAGTATTCCTTTTAATGCAAATTGGTTAAATCAACAATATGCTGAAATTATCCAAGCAATCTTATTTGATGTT

what do you want to visualize? the proteins, the GC content, some other features etc.?

Ultimately, I would like to annotate protein coding regions and compare them to other genomes.

What have you tried to do? Which code have you tried and what errors and message you got?

I would like to visualize annotated genes and compare a set of genomes.

This question seems too broad. Could you expand what do you want to compare in a set of genomes, or what type of annotation do you want to see (coding regions or proteins or ORF?)? Also please decide if your question is about seeing annotations OR compare a set of genomes (to compare genomes one would usually calculate some metrics in each genome and compare them). Remember you can [edit] the question.

It is not broad. It is very specific. I am looking for an example to plot SOMETHING with Gviz. An example, for someone who has never used it before.

I am not familiar with Gviz, but I am familiar with a variety of tools for visualizing genome features. At the scale of entire microbial genomes, the primary sequence itself is an overwhelming amount of information and can't really be visualized in a meaningful way by itself. I mean, you could color the nucleotides, but few patterns in DNA can be discerned by an unassisted human eye.
What we're typically interested in are features encoded in the genome sequence: protein-coding genes, non-coding RNAs, transposable elements, regulatory elements, binding sites, epigenetic modifications, and so on. A particular scientist is usually only interested in a subset of these features, based on what particular questions they are investigating. In your comment you say you "would like to annotate protein coding regions and compare them to other genomes", so it seems like we're on the same page here.
Visualization tools often display the genome linearly (even, for example, if it's a circular bacterial chromosome), zoomed out so that the nucleotide sequence isn't visible, and overlayed with boxes, arrows, or other glyphs representing the genes or other features of interest. Here are some examples:

This question thread contains images generated with the GenomeTools AnnotationSketch program (gt sketch command).
The UCSC genome browser which provides web-based access to dynamic displays of the human genome and many other mammalian genomes.
The IGV browser can be downloaded and easily run on your personal laptop or desktop computer.
I'm not familiar with Gviz, but it definitely looks like it belongs in this class of tools.

In each of these cases, the visualization software is not doing the work to determine the positions of the genes and other genome features. This was done independently using tools designed for this task, storing the information in an annotation file in GFF3, GTF, or BED format. The visualization tools need both the genome sequence (in Fasta format) and one or more annotation files (in GFF3, GTF, or BED format) to make meaningful visualizations.

I am touching on genomics in a project. So I am quite new to the field and don't really know what is possible, feasible. Thanks a lot for your answer. I am working with bacterial genomes.

I use prokka to annotate the genome.

Prokka should create a GFF3 file, which you can use with the Fasta file to visualize your genome. But of course, you're not the first person to annotate S. aureus, so you may also consider downloading a reference annotation for the species from NCBI RefSeq.

We will generate lots of genomes that we will have to annotate though. S aureus is not S aureus. There is a spectrum of phenotypes, they all have a different genome. And we will also not be limited to SA. The question here is not how to annotate the genome, but rather the steps necessary to visualize something.

Great. That type of background is very useful to include in the original question, so that the community knows how to focus their answers. I hope this helped!

I added it to the main question. Thanks for your comment.


text:
Ploting FDR along with the pathway as heatmap any simple way
I have this data from one of the cluster which show pathway in context of those genes but when i do a bit of filtering i see plenty of pathways with same or similar FDR which are sort of redundancy i suppose in terms of over-representation ,so i came up something like this
may be i need a bit more streamlining into it which as of now i can;t think up.
So few question just as in terms of gene clustering can i cluster those FDR or say same or similar FDR are getting clustered together which might be helpful in a way to concise and precise the result a bit better .
My data
structure(list(Function = c("covalent chromatin modification",
"histone modification", "histone methyltransferase complex",
"SWI/SNF superfamily-type complex", "histone methyltransferase activity",
"N-methyltransferase activity", "methyltransferase complex",
"protein methyltransferase activity", "histone methylation",
"S-adenosylmethionine-dependent methyltransferase activity",
"protein alkylation", "protein methylation", "macromolecule methylation",
"methylation", "methyltransferase activity", "transferase activity, transferring one-carbon groups",
"DNA helicase complex", "Ino80 complex", "peptidyl-lysine modification",
"SWI/SNF complex", "protein-DNA complex disassembly", "nucleosome disassembly",
"chromatin remodeling", "chromatin disassembly", "protein-DNA complex subunit organization",
"BAF-type complex", "INO80-type complex", "histone-lysine N-methyltransferase activity",
"histone acetyltransferase complex", "nucleosome organization",
"lysine N-methyltransferase activity", "protein-lysine N-methyltransferase activity"
), variable = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = "FDR", class = "factor"),
    value = c(1.47992665134053e-13, 1.47992665134053e-13, 7.9468967129082e-09,
    7.9468967129082e-09, 1.1841058070482e-08, 1.75996948641137e-08,
    1.75996948641137e-08, 4.71690557750932e-08, 9.63620605173175e-08,
    3.94626427029223e-07, 4.55323040178062e-07, 4.55323040178062e-07,
    4.55323040178062e-07, 8.9471201991238e-07, 1.38187416380405e-06,
    2.07900888100729e-06, 2.60629648189489e-05, 2.60629648189489e-05,
    2.91531936861022e-05, 3.19327841103331e-05, 5.04466386223372e-05,
    5.04466386223372e-05, 5.45938617344686e-05, 5.4956792922593e-05,
    5.4956792922593e-05, 5.4956792922593e-05, 8.32575963760242e-05,
    0.000287420280100368, 0.000372921199469602, 0.000387140224358431,
    0.000447883223319561, 0.000447883223319561)), row.names = c(NA,
-32L), class = "data.frame")

My code
p <- ggplot(aes(x=Function, y=variable, fill=value), data=plot.data)
fig <- p + geom_tile() + scale_fill_gradient2(low="blue", mid="white", high="red") +
  #   geom_text(aes(label=stars, color=value), size=8) + scale_colour_gradient(low="grey30", high="white", guide="none") +
  #geom_text(aes(label=stars), color="black", size=5) +
  #labs(y=NULL, x=NULL, fill="t-value") + geom_vline(xintercept=1.5, size=1.5, color="grey50") +
  theme_bw() + theme(axis.text.x=element_text(angle = -45, hjust = 0))
fig

Any help or suggestion or help would be highly appreciated

Thanks for providing the data and the code to reproduce the upper image. But how it is related to the second image below?  The redundancy you refer is that  pathways have more or less with the same genes or because they have similar FDR values? (I understand that you want to reduce from 32 pathways to less pathways)

"But how it is related to the second image below" not related perhaps that was the simpler way to convey what im trying to do . "I understand that you want to reduce from 32 pathways to less pathways)" yes that is what my goal is to ,not sure if its conceptually correct but a little search tells me there are very related may be few genes apart. Second image sorry I mean from the second image i would like to do the figure2 the middle image in the second figure where modules are plotted

The modules mean that they separated somehow the genes and still in several grouping the GO terms appear enriched. It is possible to see which pathways are similar. See this question of mine and the solutions proposed

Possible duplicate of How to select the most representative pathways from a gene enrichment analysis?

okay let me check , not sure about the duplicate but again let me go though it

looks like its not a simple visualization problem of FDR and pathway i have find out first the overlapping gene it seems there was a similar solution proposed on biostar which i had used where genes along with the pathway were plot , but this looks the second fig2 more simpler in terms of providing information overall .

Oh, you refer to the dendrogram clustering the GO terms! Yes, you can calculate the similarity of the pathways: you can use BioCor (which I am the author) for that. Afterwards order them according to their similarity

okay cool to know will try it out


text:
Changing this function to work
I have this matrix of expression data
> head(data[,1:6])
  Symbol      cl1    cl1.1     cl1.2     cl1.3    cl1.4
1   ACTB 14.07922 14.25211 13.547206 14.046814 15.33551
2 ATP5F1 10.43386 10.34219  9.671036  9.862641 10.39786
3   DDX5 13.99250 13.29734 13.706894 13.084000 13.73379
4  EEF1G 14.35303 14.26112 13.523434 13.327333 13.97747
5  GAPDH 16.07695 15.61435 15.454911 15.233670 15.90100
6    NCL 14.35149 13.32919 13.051861 12.704973 13.55680
>

I have 24 samples as cl1 and 30 samples as cl2
I want to do a t.tset, I have this function but I am not sure Why I am getting error by
volcan <- function(data, cl1, cl2) {
   tt <- t(apply(data[, -1], 1, function(d){z <- t.test(d[cl1], d[cl2]); c(z$estimate, z$p.value)}))
     tt[, 1] <- tt[, 1] - tt[, 2]
     tt <- tt[, -2]
     tt <- cbind(symbol=data[, 1], as.data.frame(tt), rank=NA)
     tt[, 3] <- -log(tt[, 3], 10)
     colnames(tt)[2:3] <- c("log.ratio", "neg.log(t.test)")
     tt[, 4] <- sqrt(tt[, 2]^4  tt[, 3]^2)
     tt <- tt[order(tt[, 4], decreasing=T), ]
     tt
}
    Error: unexpected symbol in "tt[, 4] <- sqrt(tt[, 2]^4 tt"

When I am running function as whole says that
Error: unexpected '}' in "}"

When I am putting

tt= volcan(data,cl1,cl2)

and doing line by line going down after
tt[, 4] <- sqrt(tt[, 2]^4  tt[, 3]^2)

Says
Error: unexpected symbol in "tt[, 4] <- sqrt(tt[, 2]^4  tt"

I am not sure why, I put comma in
 tt[, 4] <- sqrt(tt[, 2]^4 , tt[, 3]^2)

But another error

EDIT: I'm editing my answer as I don't have comment privileges and don't wish to add multiple answers.
Please don't change the entire premise of the question after you get an answer. If you do, retain the original content so the post makes some sense to people.
Run ?sqrt and check how the function is to be invoked, and how you're invoking it. Does it even accept 2 input parameters?

Original Content
You are subsetting the 2D object data improperly. How can data be subset as both data[,1:6] and d[data[,1:24]] where d = data? Fix that and your r code will work OK (syntactically at least).

Sorry, instead how I can define c1l and cl2 classes for first function so I don't have to subset data manually in second function

Use grepl on the colnames(d) to pick relevant columns and then subset by the results of the grepl. You should really Google this stuff.

Hang on, in the first case you defined a function with 3 parameters and used just 1. What gives?

Thank you but your comment about grep I could run the function and I edited my post

Why am I getting error - the error is pretty explicit. The function is not being invoked properly.


text:
Complexheatmap update issue
After i updated my complexheatmap library i can;t use my old annotation to label sample this is the image how i used to get something like this

Code for the below image
require(ComplexHeatmap)
require(circlize)
require(cluster)
df <- read.csv("heatmap_folder/NORMAL_EPI_TF1.txt",header = TRUE,sep = "\t")

head(df)
tail(df)
dim(df)
heat <- t(scale(t(df[,3:ncol(df)])))
head(heat)

rownames(heat) <- df$Gene
Sample = c(rep("HSC",4),rep("CMP",4),rep("GMP",4),rep("Mono",4))

anno <- as.data.frame(Sample)
categories <- data.frame(anno)
colnames(anno) <- c("Condition")

ColAnn <- HeatmapAnnotation(df=anno, which="col", col=list(Condition=c("HSC"="black","CMP"="brown",
                                                                       "GMP"="maroon1","Mono" = "darkgreen")))
RowAnn <- data.frame(df$Family)
colnames(RowAnn) <- c("Gene family")
RowAnn <- HeatmapAnnotation(df=RowAnn, col=colours, which="row")

boxAnnCol <- HeatmapAnnotation(boxplot=anno_boxplot(heat, border=TRUE, gp=gpar(fill="#CCCCCC"),
                                                    pch=".", size=unit(2, "mm"), axis=TRUE, axis_side=NULL, axis_gp=gpar(fontsize=10)),
                               annotation_width=unit(c(1, 6.5), "cm"))

boxAnnRow <- rowAnnotation(boxplot=row_anno_boxplot(heat, border=TRUE, gp=gpar(fill="#CCCCCC"),
                                                    pch=".", size=unit(3, "cm"), axis=TRUE, axis_side="top", axis_gp=gpar(fontsize=12)),
                           annotation_width=unit(c(3), "cm"))

myCol <- colorRampPalette(c("navyblue", "white", "red"))(100)
myBreaks <- seq(-2, 2, length.out=100)
hmap <- Heatmap(heat,
                name="Z-score",
                col=colorRamp2(myBreaks, myCol),
                heatmap_legend_param=list(color_bar="continuous",
                                          legend_direction="vertical", legend_width=unit(5,"cm"),
                                          title_position="topcenter", title_gp=gpar(fontsize=10, fontface="bold")),

                #Split heatmap rows by gene family
                split=df$Family,
                #split = NULL,

                #Row annotation configurations
                cluster_rows=TRUE,
                show_row_dend=FALSE,
                row_dend_side = "left",
                #row_title="Transcript", #overridden by 'split' it seems
                row_title_side="left",
                row_title_gp=gpar(fontsize=8,fontface="bold"),
                row_names_gp = gpar(fontsize = 8),
                show_row_names=TRUE,
                row_names_side="right",
                row_title_rot=0,

                #Column annotation configuratiions
                cluster_columns=FALSE,
                show_column_dend=TRUE,
                column_title="Samples",
                column_title_side="top",
                column_title_gp=gpar(fontsize=10, fontface="bold"),
                column_title_rot=0,
                column_names_gp = gpar(fontsize = 20, fontface="bold"),
                show_column_names=TRUE,

                #Dendrogram configurations: columns
                clustering_distance_columns="euclidean",
                clustering_method_columns="complete",
                column_dend_height=unit(10,"mm"),

                #Dendrogram configurations: rows
                clustering_distance_rows="euclidean",
                clustering_method_rows="complete",
                row_dend_width=unit(15,"mm"),

                top_annotation_height=unit(0.5,"cm"),
                top_annotation=ColAnn,

                bottom_annotation_height=unit(2, "cm"),
                bottom_annotation=boxAnnCol)
hmap + RowAnn

The code remains the except the data
Now not sure may be im messing up the annotation so that i can;t see sample label and the family annotation
Since this was the error

"Error: axis_gp is removed from the arguments. Use axis_param =
  list(gp = ...) instead." which i modified a bit

Then i modified my code to this
RowAnn <- data.frame(df$Family)
colnames(RowAnn) <- c("Gene family")
RowAnn <- HeatmapAnnotation(df=RowAnn, col=colours, which="row")

ColAnn <- HeatmapAnnotation(df=anno, which="col", col=list(Condition=c("HSC"="black","CMP"="brown",
                                                                       "GMP"="maroon1","Mono" = "darkgreen")))

boxplotCol <- HeatmapAnnotation(boxplot=anno_boxplot(heat,
  border=TRUE, gp=gpar(fill="#CCCCCC"), pch=".",
  size=unit(2, "mm"), axis=TRUE, axis_param = list(side = "left") ,
  annotation_width=unit(c(1, 5.0), "cm"), which="col")
)

hmap <- Heatmap(heat, name = "expression",col=colorRamp2(myBreaks, myCol),
                split=df$Family,
                #split = NULL,

                show_row_names = TRUE, show_column_names = TRUE, cluster_rows = TRUE,
                cluster_columns = TRUE, show_column_dend = TRUE, show_row_dend = TRUE,
                row_dend_reorder = TRUE, column_dend_reorder = TRUE, clustering_method_rows = "ward.D2",
                clustering_method_columns = "ward.D2", width = unit(100, "mm"),

                #top_annotation_height=unit(1.0,"cm"), top_annotation=colAnn,

                bottom_annotation=boxplotCol)

draw(hmap+RowAnn, heatmap_legend_side="left", annotation_legend_side="right")

my fig
I would be glad to know  what im doing wrong as i can;t she condition getting annotated
my data
dput(df)
structure(list(Family = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("EPI",
"TF"), class = "factor"), Gene = structure(c(31L, 23L, 32L, 40L,
26L, 34L, 14L, 17L, 39L, 21L, 15L, 30L, 9L, 1L, 33L, 2L, 5L,
11L, 18L, 16L, 7L, 27L, 8L, 19L, 22L, 12L, 24L, 36L, 6L, 20L,
28L, 41L, 43L, 44L, 45L, 47L, 48L, 49L, 50L, 13L, 17L, 38L, 42L,
21L, 25L, 35L, 37L, 10L, 11L, 4L, 46L, 3L, 29L), .Label = c("BABAM1",
"BAP1", "BRF2", "CASZ1", "CBX8", "CCDC17", "CEBPA", "CUX2", "DAPK3",
"DEAF1", "DOT1L", "EBF1", "ERF", "HCFC1", "HDAC7", "HES6", "HMG20B",
"INO80E", "IRF8", "KLF16", "MBD3", "MEF2D", "MTA3", "MXD4", "PIN1",
"PIWIL4", "PREB", "REPIN1", "RFX2", "SETD1A", "SFMBT1", "SFMBT2",
"SMARCA4", "SMYD2", "SNAPC4", "SREBF2", "TCF20", "TCF7", "TFPT",
"UBE2A", "ZBTB17", "ZFHX2", "ZNF142", "ZNF319", "ZNF414", "ZNF446",
"ZNF469", "ZNF581", "ZNF646", "ZNF76"), class = "factor"), HSC1 = c(1.38256100990102,
6.03508433561645, 8.53718694596963, 9.7579, 4.9089643439748,
6.08652168925513, 10.85362, 5.791612, 4.717658, 5.034092, 7.6573073927514,
8.214206, 6.087519, 7.795752, 9.147491, 9.16065391580272, 5.475244,
6.839048, 4.307532, 0.4008883181761, 3.681525, 5.39357937244362,
0.048702874829245, 4.884674, 9.42704096616912, 0.331449426359954,
8.70131868583511, 9.860997, 5.96926894038996, 5.122815, 9.84364182958174,
4.99646644487682, 8.64613509396529, 6.221252, 4.788327, 4.221885,
6.904238, 7.027222, 6.90290115677703, 5.910114, 5.791612, 7.85232611449929,
4.30633834185685, 5.034092, 6.587824, 3.77111448331312, 4.816565,
5.83621984608683, 6.839048, 5.38002542796833, 2.499083, 4.07296733289583,
7.384605), HSC2 = c(3.3250440950498, 5.69710569315867, 7.53535864842004,
9.46294, 4.63015366085456, 6.27268560573257, 10.61808, 5.936903,
5.754081, 5.17955, 7.34546746135814, 8.156995, 4.959736, 7.1724,
9.519366, 8.399964626033, 3.79417, 6.580323, 3.006501, 0.465543296427136,
3.854056, 6.19801925234426, 1.35661238539194, 2.520579, 10.1179535946398,
0.438637606305892, 8.02299949380521, 9.599955, 5.14592215971399,
5.362818, 9.62067506274524, 5.81177950264385, 8.6388082408386,
6.580456, 5.507812, 4.647564, 5.933962, 7.036506, 6.61438422233323,
7.081847, 5.936903, 7.25061985280206, 6.14119015507351, 5.17955,
5.51992, 4.7980903060972, 4.688224, 5.75480261400296, 6.580323,
4.61480421865767, 2.473632, 2.58901029457056, 7.141953), HSC3 = c(4.5609308176212,
6.31335736623603, 8.08532370258223, 9.954938, 5.33971874563167,
6.31519305036765, 8.686445, 3.652774, 1.705097, 3.590782, 5.82598582988908,
7.160368, 2.730684, 6.544555, 8.308113, 7.28529068473834, 2.181592,
6.153712, 3.769699, 0.686412239668329, 3.181293, 5.81660138585872,
0.513314099277583, 3.807037, 9.40253894125251, 0.789263521770297,
6.25643885795079, 8.4667, 3.59885151800005, 1.959458, 8.16946518123644,
3.78892554507724, 8.17931637873862, 6.871673, 2.207662, 2.586419,
4.611411, 6.342816, 6.24601538582422, 4.504617, 3.652774, 6.1448035804027,
4.32474896196878, 3.590782, 4.283432, 3.97377694799887, 5.081613,
4.67353303760197, 6.153712, 2.95856351440761, 1.797669, 2.61144324906345,
6.927889), HSC4 = c(4.33336593087016, 5.51218211320577, 7.98927439884575,
10.02818, 4.72376095157983, 7.06408657043682, 10.6532, 6.332262,
4.356266, 5.063746, 7.40605558817077, 7.943625, 4.577773, 7.114057,
10.10959, 8.39439253099313, 0.843862, 7.275198, 3.125764, 0.461081947661647,
4.914624, 6.47822557945205, 0.151047922603591, 4.233378, 9.76778876120242,
1.52153581933298, 7.74186897540428, 9.848883, 4.27034862337376,
5.461201, 9.20573611283668, 5.57160015049504, 8.90499660789483,
5.972766, 4.37926, 4.914465, 6.288035, 6.667292, 6.48347655570048,
6.525983, 6.332262, 6.67910931337338, 5.93921207229708, 5.063746,
5.391649, 5.14382963917955, 3.682687, 6.17864627473077, 7.275198,
4.55820389425563, 2.464959, 3.36552973589941, 7.155267), CMP1 = c(6.87345497292835,
5.51412758368083, 7.31352446967914, 10.55862, 6.25444074607741,
7.17818981714724, 11.09789, 6.353011, 4.176257, 5.988031, 8.45885235860671,
8.239023, 4.831419, 8.439534, 9.741489, 9.12973766558991, 1.63498,
7.616455, 4.40003, 3.35836219130892, 6.782445, 8.09320436383016,
0.034087182831375, 7.328077, 10.1348328626404, 3.77197918672699,
8.22105235434344, 9.586611, 3.96736165755976, 6.120908, 9.66405823977144,
5.5410623284247, 7.98439042490943, 5.617969, 4.389489, 3.320773,
7.259381, 6.342179, 7.80157963446663, 6.00359, 6.353011, 4.88520066048698,
4.57336386910792, 5.988031, 7.534274, 4.22034258843835, 4.407149,
6.04554630191099, 7.616455, 5.5286242078284, 2.938658, 6.05835585199437,
7.254343), CMP2 = c(4.57632794651802, 3.78653056181795, 7.38323687896953,
10.18887, 5.42618113653324, 7.14753074148422, 11.10573, 6.444083,
5.364945, 5.56755, 7.87006213457315, 8.183806, 5.547153, 7.408484,
10.41511, 8.5734105703595, 3.332623, 7.927959, 3.044397, 2.23975823913519,
5.948004, 6.8523382099618, -0.208208601046198, 2.502613, 9.61958525630466,
0.083082776643226, 7.96816766770755, 9.658819, 4.28029266262047,
6.064432, 9.51407062066367, 6.0002887227803, 8.84860518364821,
6.63929, 5.418921, 4.246354, 6.265858, 6.726515, 6.59033513636581,
6.599091, 6.444083, 6.42468503826279, 5.39157032731654, 5.56755,
5.975197, 5.49143188875569, 4.882971, 6.3269787223519, 7.927959,
4.52014312274917, 3.02202, 3.32107044636493, 7.234948), CMP3 = c(5.85154481845215,
4.70657054545281, 7.81269743162002, 10.51041, 6.02525567584962,
7.6061510389796, 8.938033, 4.185122, 3.453336, 4.250061, 6.95084523398581,
7.016165, 3.711631, 6.691113, 9.291314, 6.79870749458327, 2.556527,
6.640901, 2.472737, 0.713500223112626, 4.610029, 7.09501723097758,
0.554991742557483, 3.663828, 8.64957961123419, 0.830920666997163,
5.03120897448563, 8.193863, 3.66385830275847, 2.514419, 8.57920717695102,
3.68106350062846, 7.99067978896729, 4.917504, 1.705394, 2.940584,
4.935565, 4.679952, 5.45403468978382, 4.097788, 4.185122, 4.7017134210513,
3.36996680968093, 4.250061, 4.841303, 4.16528140537129, 4.311709,
4.34280333979994, 6.640901, 3.46696931369128, 1.836416, 1.2954688478789,
6.759544), CMP4 = c(5.70607301130076, 4.90995364796149, 8.07668148492955,
10.49127, 6.12272998130959, 8.1070578791571, 11.14272, 6.069435,
3.77716, 5.905144, 7.80284615350651, 8.193236, 5.266133, 7.994962,
10.62971, 8.76622729912917, 1.728943, 8.018743, 2.64787, 3.17362412630739,
6.675791, 8.17425862958478, 1.3062705649262, 8.305058, 9.60100662458438,
3.64083547966437, 7.03151917955059, 9.697494, 4.22387311586196,
5.872986, 9.48040032459999, 5.9054376735385, 9.0993222961645,
6.641411, 4.042785, 4.175351, 6.040673, 7.025595, 6.8701243794046,
6.474802, 6.069435, 5.98875791029547, 5.44536133412665, 5.905144,
6.253757, 5.55713735563886, 3.183721, 6.23569788279688, 8.018743,
4.14168543985438, 3.140228, 2.85944735485205, 7.19532), GMP1 = c(6.00003380096822,
6.45782274084263, 10.1386873750581, 11.18212, 4.99574155426759,
8.22548270576452, 11.80921, 6.139842, 4.716938, 6.208399, 7.25702470622416,
9.216353, 4.459075, 7.604192, 10.214, 9.71979154806924, 4.045693,
6.892321, 4.36733, 4.54209088244969, 4.363115, 7.42142574723769,
8.16658165526266, 13.13801, 11.7030896033115, 1.38480651276548,
7.70549837966537, 10.51655, 5.6799276696584, 7.117451, 10.4613436109924,
7.51060888343217, 9.45360701936786, 11.02082, 4.741664, 4.37359,
6.757818, 7.101331, 6.51838172412493, 6.234017, 6.139842, 7.16147421243415,
5.80674368005962, 6.208399, 6.860234, 4.69772329684801, 5.404108,
6.45675195710298, 6.892321, 5.30230900977469, 4.155601, 4.36874064027233,
7.42003), GMP2 = c(5.74747276487971, 5.97553824097351, 8.84397150139003,
10.54739, 7.01942050935297, 8.32418895763336, 11.28484, 6.100887,
5.371866, 5.782092, 6.81459497341517, 8.300889, 4.573727, 7.155416,
10.55908, 8.53584780642442, 3.288921, 7.981419, 3.343707, 1.78160205563733,
7.106303, 6.97534663269409, 4.30909405726115, 10.94426, 10.9546928815164,
6.74862372900741, 7.1870575783895, 9.650893, 4.16099729988147,
6.067064, 9.97780340347263, 5.92239739476255, 9.00840651100957,
5.952921, 5.034178, 4.255218, 5.533772, 6.671744, 6.58968157482356,
6.572507, 6.100887, 6.06061888187195, 4.67987580614607, 5.782092,
5.976481, 5.57605980178547, 4.822138, 5.9343704738438, 7.981419,
4.0036610231274, 2.709241, 2.35084397835661, 7.122418), GMP3 = c(6.86258223460792,
6.20078322375291, 9.61446552902509, 10.80572, 7.36372917472567,
8.3826076358515, 11.15008, 5.201591, 4.123945, 5.045783, 6.24154676477513,
8.276453, 4.263376, 7.519041, 10.51785, 8.30200623104098, 2.495567,
7.857207, 2.962492, 1.93915883789869, 7.249551, 7.84810925448229,
5.84602007858824, 11.77414, 11.1027524778454, 5.63020389504886,
6.50727637405781, 9.456897, 3.75134646474614, 5.932853, 10.1039559568546,
4.9400044296399, 9.31472753734898, 6.039605, 3.336343, 3.975937,
5.736387, 6.905666, 6.30344093479224, 5.426484, 5.201591, 5.47026133533826,
4.69976751703472, 5.045783, 5.749752, 5.23937658540781, 5.154187,
6.34478978856221, 7.857207, 3.50805939441741, 2.657933, 2.78283797334871,
7.091083), GMP4 = c(6.05880869638876, 6.28654396576068, 9.96479543345143,
10.73202, 7.213420739218, 8.81669595423024, 11.31576, 5.679182,
4.155314, 5.673872, 6.43309255133785, 8.167664, 4.513147, 7.768381,
10.68025, 8.7651559260891, 2.485242, 8.41621, 3.457858, 2.07304058721841,
6.412872, 8.09668562132081, 5.97986866984983, 12.08623, 11.28032878957,
6.57392060927003, 6.41585754599334, 9.749641, 3.66579392385564,
5.878271, 9.75282317403084, 5.77880844881496, 9.31419691100924,
5.048616, 3.784422, 4.082659, 6.015668, 6.789027, 6.44947035574528,
6.177861, 5.679182, 5.40028170642952, 4.40100191946731, 5.673872,
6.189491, 5.32217020480945, 3.106884, 5.98621706786421, 8.41621,
3.507010215991, 2.810708, 2.26597487448895, 6.822586), Mono1 = c(5.2563684147039,
6.45633325704225, 8.35683375157621, 10.4683, 6.89578871222261,
6.89453589831897, 10.76922, 5.79888, 3.938016, 4.625604, 8.14811443428627,
8.89819, 4.565085, 7.582637, 9.690787, 9.26841770880937, 2.680727,
7.256856, 5.506176, 1.73004263683502, 7.168846, 7.34133258574044,
2.04394933289888, 10.88391, 10.8831286687177, 2.76098817217713,
7.62084289021189, 9.376334, 5.89374731019185, 5.059797, 9.58406921227736,
6.10083260496288, 8.90601784817498, 8.926934, 4.514644, 4.671412,
6.824876, 7.160325, 7.12646156474527, 6.029051, 5.79888, 7.45127447114964,
6.58727075697179, 4.625604, 6.726376, 4.54492833429885, 5.482923,
5.75789405204644, 7.256856, 5.01661617666884, 3.912327, 4.21481739920369,
7.311136), Mono2 = c(5.85774218629914, 5.24626015596457, 7.87393184409552,
10.74763, 6.1115477644807, 8.59285304183197, 11.30387, 6.487999,
5.641296, 5.605635, 6.83338037544917, 8.158746, 5.031193, 7.791207,
10.12365, 8.83146107303257, 3.839879, 7.995034, 3.864015, 1.95600811707846,
6.547528, 7.15767104450585, 1.13087340949285, 10.60541, 10.6662569861529,
2.45950067788357, 7.80075639569934, 9.772267, 4.10768104209722,
6.039518, 9.45126490091168, 5.43166366925388, 8.53913541490463,
6.852944, 5.210281, 3.909477, 6.10825, 6.759114, 6.10127101835935,
6.394382, 6.487999, 6.77443082143061, 4.65095032762504, 5.605635,
6.292733, 5.06179443435785, 4.687693, 5.97873952456647, 7.995034,
4.53685791636727, 2.036068, 3.79175025763666, 7.80933), Mono3 = c(5.79480048210947,
5.93841021449181, 8.4456428811783, 10.96997, 5.90368127873194,
7.57633522462406, 10.77521, 5.357294, 3.631277, 5.573695, 6.02937219504517,
7.823033, 4.89389, 7.134801, 9.942904, 8.29297682562671, 2.376142,
7.942619, 2.938838, 2.58663314475541, 5.473615, 6.54736080917345,
1.50818476189011, 10.39009, 9.81896143227971, 2.61179572528994,
7.11294273571537, 9.321912, 4.59467076052555, 4.846162, 9.0889984655259,
5.5580711224014, 8.60684080318993, 6.807437, 3.781213, 3.349656,
5.352098, 6.732808, 6.18873397274794, 5.580769, 5.357294, 6.02896733254353,
4.09719798968178, 5.573695, 5.576249, 4.76264681239467, 4.599397,
5.85107919870538, 7.942619, 4.46991267319193, 2.359072, 2.49361466378781,
7.168415), Mono4 = c(5.79136484133828, 5.22885638262284, 8.66188302349651,
10.80097, 6.02816710441609, 8.28180453176157, 10.95579, 5.606963,
4.781397, 6.146859, 7.35941757150855, 8.279649, 5.035288, 8.131193,
10.46082, 8.7442177262891, 2.904981, 8.009644, 2.855958, 2.31248706248342,
7.171729, 7.85215015030227, 1.5674747074644, 10.34766, 10.4041048711624,
1.7198811173707, 7.7219124608883, 9.652986, 4.80881131147327,
5.631832, 9.37540653744023, 5.93520028734304, 8.83460813839155,
7.145424, 4.609511, 4.550807, 6.666205, 6.956604, 6.66113029874611,
6.54366, 5.606963, 5.47022342513485, 4.59138757680696, 6.146859,
5.969648, 5.20105550006015, 3.857023, 6.22204033269913, 8.009644,
4.44000233834592, 2.854773, 2.88863557103663, 6.788677)), class = "data.frame", row.names = c(NA,
-53L))

Did you read the news of the package to know if anything else has changed? Did you look at the help pages to see how is the new way to label the heatmap?

I updated to use it for oncoprint after that i see some of the old codes are not working , im slowly making changes while i read it ...

I see that you created the ColAnn variable with the Condition annotation but did not include it in the draw call (or even in the Heatmap() function). Could that be the problem here?
EDIT: Oops, just noticed that this post is really old. Doesn't have an answer, so I guess I'll let mine stay.

yeah that was the mistake i thought after updating ,it broke


text:
How to analyze and visualize CRISPR-CAS9 screening results?
I have a whole-gene CRISPR screening of an assay at different time points. I created a table of the gene number and the count of its reads:
gene_id, count_t0, count_t1, count_t3
ABCA1       3        6          7

What would be a standard way to process such type of data to inference information on up-regulated/down-regulated genes.
I saw this tool: http://crispr-analyzer.dkfz.de/, but it requires my NGS raw data and these file are very large in my case. I'd expect that the information for analysis purposes can be derived solely from the counting table above.

Have you tried DESeq2? It's one of the options from the tool you linked to even.


text:
'Wildcards' object has no attribute 'sample'
I'm using Snakemake and I'm trying to get all my QC reports into multiqc but I get the following error:
WorkflowError in line 120 of /rst1/2017-0205_illuminaseq/scratch/swo-406/test_snakemake_full/Snakefile:
'Wildcards' object has no attribute 'sample'

Rule 120 is the start of my multiqc rule which looks like this:
rule multiqc:
    input:
        expand([os.path.join(analyzed_dir, '{sample}.stat', '{sample}.cnt'),
                os.path.join(rseqc_dir, '{sample}.bam_stat.txt'),
                os.path.join(rseqc_dir, '{sample}.clipping_profile.xls'),
                os.path.join(rseqc_dir, '{sample}.deletion_profile.txt'),
                os.path.join(rseqc_dir, '{sample}.infer_experiment.txt'),
                os.path.join(rseqc_dir, '{sample}.geneBodyCoverage.txt'),
                os.path.join(rseqc_dir, '{sample}.inner_distance.txt'),
                os.path.join(rseqc_dir, '{sample}.insertion_profile.xls'),
                os.path.join(rseqc_dir, '{sample}.junction.xls'),
                os.path.join(rseqc_dir, '{sample}.junctionSaturation_plot.r'),
                os.path.join(rseqc_dir, '{sample}.mismatch_profile.xls'),
                os.path.join(rseqc_dir, '{sample}.read_distribution.txt'),
                os.path.join(rseqc_dir, '{sample}.pos.DupRate.xls'),
                os.path.join(rseqc_dir, '{sample}.seq.DupRate.xls'),
                os.path.join(rseqc_dir, '{sample}.GC.xls'),
                os.path.join(rseqc_dir, '{sample}.NVC.xls'),
                os.path.join(rseqc_dir, '{sample}.qual.r'),
                os.path.join(rseqc_dir, '{sample}.RNA_fragment_size.txt'),
                os.path.join(rseqc_dir, '{sample}.STAR.genome.sorted.summary.txt'),
                os.path.join(fastqc_dir, '{sample}_R{read}_001_fastqc.html')],
                sample=samples['samples'], read=reads)
    output:
        os.path.join(qc_dir, 'multiqc_report.html')
    conda:
        "envs/multiqc.yaml"
    shell:
        '''
        #!/bin/bash
        multiqc . --outdir qc --ignore .snakemake --force
        '''

My "rule all looks like this:"
rule all:
    input:
        expand([os.path.join(analyzed_dir, '{sample}.genes.results'),
                os.path.join(qc_dir, 'multiqc_report.html')],
                sample=samples['samples'])

I don't see anything wrong with it, in fact, it is very similar to my "rule all" I used before and which worked. Can anyone see an error?
EDIT
I limited the rule to only one input, I still have the same error:
rule multiqc:
    input:
        expand(os.path.join(rseqc_dir, '{sample}.geneBodyCoverage.txt'),
                sample=samples['samples'])
    output:
        os.path.join(qc_dir, 'multiqc_report.html')
    conda:
        "envs/multiqc.yaml"
    shell:
        '''
        #!/bin/bash
        multiqc . --outdir qc --ignore .snakemake --force
        '''

This is the ouput, you can see it read the samples correctly:
run_snakemake.sh
/rst1/2017-0205_illuminaseq/scratch/swo-406/snakemake/lib/python3.6/site-packages/snakemake/workflow.py:12: FutureWarning: read_table is deprecated, use read_csv instead, passing sep='\t'.
  from collections import OrderedDict
Building DAG of jobs...
/rst1/2017-0205_illuminaseq/scratch/swo-406/snakemake/lib/python3.6/site-packages/snakemake/workflow.py:12: FutureWarning: read_table is deprecated, use read_csv instead, passing sep='\t'.
  from collections import OrderedDict
Building DAG of jobs...
Using shell: /bin/bash
Provided cluster nodes: 256
Job counts:
    count   jobs
    1   all
    1   multiqc
    2

[Mon Jan 14 17:46:40 2019]
rule multiqc:
    input: qc/rseqc/0067_P2018SEQE27R01_S1.geneBodyCoverage.txt, qc/rseqc/0067_P2018SEQE27R02_S2.geneBodyCoverage.txt, qc/rseqc/0067_P2018SEQE27R03_S3.geneBodyCoverage.txt, qc/rseqc/0067_P2018SEQE27R04_S4.geneBodyCoverage.txt, qc/rseqc/0067_P2018SEQE27R05_S5.geneBodyCoverage.txt, qc/rseqc/0067_P2018SEQE27R06_S6.geneBodyCoverage.txt, qc/rseqc/0067_P2018SEQE27R07_S7.geneBodyCoverage.txt, qc/rseqc/0067_P2018SEQE27R08_S8.geneBodyCoverage.txt
    output: qc/multiqc_report.html
    jobid: 2

WorkflowError in line 119 of /rst1/2017-0205_illuminaseq/scratch/swo-406/test_snakemake_full/Snakefile:
'Wildcards' object has no attribute 'sample'
17:46 nlv24077@kiato /rst1/2017-0205_illuminaseq/scratch/swo-406/test_snakemake_full > ll qc/rseqc/0067_P2018SEQE27R07_S7.geneBodyCoverage.txt
-rw-rw-r-- 1 nlv24077 2017-0205_illuminaseq 303 Jan  7 16:00 qc/rseqc/0067_P2018SEQE27R07_S7.geneBodyCoverage.txt

EDIT 2:
This is how I start Snakemake:
> cat run_snakemake.sh
source activate /rst1/2017-0205_illuminaseq/scratch/swo-406/snakemake
snakemake --dag | dot -T svg > dag.svg
snakemake --cluster-config cluster.json --cluster "qsub -l nodes={cluster.nodes}:ppn={cluster.ppn} -N {wildcards.sample}_{rule}" --jobs 256 --use-conda --rerun-incomplete

EDIT 3:
Ok, that last remark (by @finswimmer) triggered me, I removed the "{wildcards.sample}_" from the above command and then it works. Or, works, I get a very different error, it is multiqc specific and solved as specified here: https://multiqc.info/docs/#locale-error-messages. When adding the two lines as suggested this error disappears and everything works as expected, even with the long inputs list.
RuntimeError: Click will abort further execution because Python 3 was configured to use ASCII as encoding for the environment. Consult https://click.palletsprojects.com/en/7.x/python3/ for mitigation steps.

Anyway, I'm guessing Snakemake does not understand how to name this task when it uses an entire list of {sample}'s... Is this a bug? What should be used? I guess it is impossible to use {sample} for naming for this rule. It's a shame I have to kill the naming for all rules because using sample_rule is very nice for keeping track.

Can you start by limiting the input to that rule to a couple files? That will at least help debug this.

I don't know anything about snakemake but the problem seems to come from not declaring what a sample should be

@llrs sample is defined in the expand() function. Oddly, there's no direct use of wildcards outside of that, which is usually where this sort of thing crops up.

Is samples['samples'] defined somewhere?

And how does your command look like to start the workflow?

I think you could use it, but not that way. But I think that would be another and please post how you solved it (I want to learn to use snakemake).

Hi @llrs, the answer is in EDIT 2, there is my "start snakemake" script. I used the wildcards.sample to name the jobs on our qsub cluster, but I can't do that because my MultiQC does not have a single defined wildcards.sample (it has many). So I simply removed that part (the {wildcards.sample}) and now I only have the rules as the jobname, sadly. Does that help? How do you think I could use it then?

@Freek Well, to signal that you got an answer here (on the stackexchange.com sites) we post an answer! You can post an answer to your own question (it is completely normal and encouraged). In the answer please explain how it solved your problem as you did in this comment.

Ok, thanx, I'll remember that, for now James Hawley worded it nicely already :)

Snakemake gets wildcards from the parsing the input/output file names.
For example, if you had
rule sort:
    input: "{file}.bed"
    output: "{file}.sorted.bed"
    shell: "sort -k1,1 -k2,2n {input} > {output}"

then the wildcards variable would be set as wildcards = {"file": ...}.
In your multiqc rule, you have the {sample} variable within your expand(), but this variable doesn't exist outside of expand's local context. You don't actually have any wildcards there, nor in your output file (qc_dir and rseqc_dir are variables you define elsewhere, but not "wildcards" as Snakemake interprets them).
So you're getting the error 'Wildcards' object has no attribute 'sample' because there are no wildcards for that rule.
Your rule is currently set up to process every sample that you have listed in  samples['sample']. If you want your multiqc rule to work for a single sample, simply remove the sample=samples['sample'] and double brace the {sample} wildcard.
rule multiqc:
    input:
        expand([os.path.join(analyzed_dir, '{{sample}}.stat', '{{sample}}.cnt'),
                os.path.join(rseqc_dir, '{{sample}}.bam_stat.txt'),
                os.path.join(rseqc_dir, '{{sample}}.clipping_profile.xls'),
                os.path.join(rseqc_dir, '{{sample}}.deletion_profile.txt'),
                os.path.join(rseqc_dir, '{{sample}}.infer_experiment.txt'),
                os.path.join(rseqc_dir, '{{sample}}.geneBodyCoverage.txt'),
                os.path.join(rseqc_dir, '{{sample}}.inner_distance.txt'),
                os.path.join(rseqc_dir, '{{sample}}.insertion_profile.xls'),
                os.path.join(rseqc_dir, '{{sample}}.junction.xls'),
                os.path.join(rseqc_dir, '{{sample}}.junctionSaturation_plot.r'),
                os.path.join(rseqc_dir, '{{sample}}.mismatch_profile.xls'),
                os.path.join(rseqc_dir, '{{sample}}.read_distribution.txt'),
                os.path.join(rseqc_dir, '{{sample}}.pos.DupRate.xls'),
                os.path.join(rseqc_dir, '{{sample}}.seq.DupRate.xls'),
                os.path.join(rseqc_dir, '{{sample}}.GC.xls'),
                os.path.join(rseqc_dir, '{{sample}}.NVC.xls'),
                os.path.join(rseqc_dir, '{{sample}}.qual.r'),
                os.path.join(rseqc_dir, '{{sample}}.RNA_fragment_size.txt'),
                os.path.join(rseqc_dir, '{{sample}}.STAR.genome.sorted.summary.txt'),
                os.path.join(fastqc_dir, '{{sample}}_R{read}_001_fastqc.html')],
                read=reads)
    output:
        # this {sample} is in single {} and is recognized as a wildcard
        os.path.join(qc_dir, '{sample}', 'multiqc_report.html')
    conda:
        "envs/multiqc.yaml"
    shell:
        '''
        #!/bin/bash
        # using {wildcards.sample} here since MultiQC likes directories
        multiqc . --outdir qc/{wildcards.sample} --ignore .snakemake --force
        '''

When wildcards are recognized, you'll see them in the job list.
Running snakemake 0067_P2018SEQE27R01_S1/multiqc_report.html, for example, should give you something like:
rule multiqc:
    input: qc/rseqc/0067_P2018SEQE27R01_S1.geneBodyCoverage.txt, ...
    output: qc/0067_P2018SEQE27R01_S1/multiqc_report.html
    jobid: 2
    wildcards: sample=0067_P2018SEQE27R01_S1

Thank you, nicely explained, I'll accept this answer. But, do you know of any way to temporarily set the wildcards.sample variable? So that I can use it for naming my multiqc job? I find it a shame that because of this one rule I can't have the {sample}_{rule} naming scheme anymore...

Maybe you could try using multiple bash commands? Something like multiqc . --outdir qc --ignore .snakemake --force && mv multiqc_report.html {wildcards.sample}_multiqc.html?

But I think that because multiqc scans the entire directory, all your files are going to get scanned anyway, and end up in each report. So I don't know if they're actually going to be any different, or if you're just going to end up with the same file under multiple names

@JamesHawley, I thought the OP wanted many inputs (defined by the expand() call) to a single output (qc/multiqc_report.html), rather than your solution, which if I understand results in an output for every input. It should be possible to achieve what the OP asked for, and the multiqc wrapper seems to suggest so. In which case, what was wrong with the OP's code? (I'm having the same error in a similar situation, thus my comment)

I got it now: the problem wasn't with the rule's code itself, it was to do with the call to --cluster. So, there is no known solution for this problem of wanting to use the wildcards.sample in the call to qsub when there is one, but otherwise not? This would be particularly useful for log files.

Although @JamesHawley answer was informative, I don't think it solved the initial issue requested by the OP. (edit: or maybe it partially did, I'm not sure as I got a bit confused with the wildcards on the initial question/answer. I'm anyway leaving this answer here, in case another way of explaining it is useful for others)
Solution
I think you only need to change the output of your rule multiqc from:
os.path.join(qc_dir, 'multiqc_report.html')
to
os.path.join(qc_dir, '{sample}.html')
This seemed counter-intuitive to me initially, but I think what's happening is that snakemake will look at rule all and see that you want an output called "qc/multiqc_report.html". And the only place where it can find something of that sort is in rule multiqc, by replacing the wildcard {sample} by the string multiqc_report.
It felt confusing to me, because the same {sample} wildcard is used in other rules (namely the rules you use to generate the outputs in your analysed_dir). But snakemake will be smart to replace {sample} with whatever it needs to recreate the output files you're requesting.

Longer explanation
Here is a (hopefully) reproducible example, if it helps.
First create some fake directory structure with some files:
mkdir inputs
mkdir analysed_dir
mkdir qc_dir

touch inputs/file1.txt
touch inputs/file2.txt

Then, this would be the snakefile:
# Define your variables, fetch things from config, or csv, etc...
input_dir = "inputs"
output_dir = "outputs"
samples = ["file1.txt", "file2.txt"]

# Rule all defining all the output files required
rule all:
  input:
    expand("analysed_dir/{samples}", samples = samples),
    "qc_dir/multi_qc_report.html"

# this rule is equivalent to whichever rules process your data individually
rule analyse_each_file:
  input:
    "inputs/{sample}.txt"
  output:
    "analysed_dir/{sample}.txt"
  shell:
    "cat {input} > {output}"

# This rule is equivalent to your multiqc step where multiple files result in a single output
rule multiqc:
  input:
    expand("analysed_dir/{samples}", samples = samples)
  output:
    "qc_dir/{sample}.html"
  shell:
    "cat {input} > {output}"

to generate this pipeline:

And then issuing:
snakemake --cluster "qsub -N {wildcards.sample}-{rule}" --jobs 10
Generates the following job names:

file1-analyse_each_file
file2-analyse_each_file
multi_qc_report-multiqc

Edit with further explanation:
Notice in the DAG how the sample wildcard takes different values for the different rules. I think the snakemake logic is as follows:

Starting with rule all I see that you want three files as input (the first two come from the expand() call):

analysed_dir/file1.txt
analysed_dir/file2.txt
qc_dir/multi_qc_report.html

Let me start with the first one analysed_dir/file1.txt. The file doesn't exist, so let me look in the other rules and see if I can find it in their output. Actually in rule analyse_each_file I see the output is set to "analysed_dir/{sample}.txt", so if I replace {sample} by file1.txt I can recreate your file from it. Therefore, the input for that rule will be "inputs/file1.txt" (and this file already exists, so I can issue this job).
For the second file analysed_dir/file2.txt, well it's a similar case as with the previous one, so I'll issue a job with the same rule analyse_each_file, but this time {sample} = file2.txt.
Finally, you want a file called qc_dir/multi_qc_report.html. Where can I create such a file from? Well, the output of rule multiqc is set to "qc_dir/{sample}.html", so if I make {sample} = multi_qc_report, then I can issue the job and produce your file of interest. In this case, the rule takes as input two files "analysed_dir/file1.txt" and "analysed_dir/file2.txt", which I already know how to create from another rule.

The key thing to realise is that snakemake defines the wildcard {sample} for each rule separately. So, in this example, wildcard.sample in rule analyse_each_file does not take on the same values as in rule multiqc.

Thank you for your comment, I still don't have a satisfying solution (meaning I now just leave out the sample wildcard from the qsub job names). For your solution I wonder, will Snakemake not complain about "qc_dir/{sample}.html" never being generated? I think the core issue is that the Multiqc rule takes in all samples so it is correct that there is no single "sample" wildcard. I was just looking for a way to set the sample wildcards to something like "all_samples" for the purpose of naming the qsubbed jobs. Am I making any sense?

@Freek, when you say "I still don't have a satisfying solution", do you mean you tried my solution and it didn't work? (note, I've edited my solution to more correctly match your posted code using os.path.join() - see the Solution section at the top). I don't think snakemake will get confused, because it seems to me that wildcards are defined within each rule individually. See my further edit with expanded explanation at the end of the answer (below the DAG). I hope it helps.


text:
log rank p value interpretation from gepia 2
How do i interpret the overall survival when the p value reported is not significant going by the standard convention .
Not sure what kind of test it runs while comparing the expression may be one non-parametric test , but the survival is kaplan meier on cox ph model which they have mentioned .
So is it a significant one if i have to tell in context of survival ?

Could you provide a bit of context? From which article did you get this? Why do you think it needs different interpretation when the p-value is above your threshold?

it's from gepia site so im looking for genes in databases which are significant across patient sample ,boxplot shows significant but the survival is not significant if i understand ,why is it so?

The plot on the left shows that the gene has very different expression between two groups. However, according to the graph on the right, these high expression is not directly related to survival because both groups have similar survival time.
The gene expression so high might be to compensate for other deficiencies or problems.

yes both of them have similar survival  which is reflected in the log rank value ? in other words if i see the log rank value is not significant ,then i shouldn't consider it..

@krushnach You shouldn't consider it for what? Please [edit] the question to include more details (I don't even know which gene is: LAML or zbtb16??), perhaps a link to the result on the website is possible

okay i will do it. i cannot link it as it the result is stored and rendered there but here is the link http://gepia2.cancer-pku.cn/#survival the data set im using LAML gene is ZBTB16 you can have a look..

Ok, that's for the survival curve and the other plot were does it come from?

http://gepia2.cancer-pku.cn/#analysis this one

Ok, know I know how did you get these plots. But why do you need them ? Why are you looking at these plots? What are you looking for

as they were getting differentially expressed in leukemic stem cells and blast cells so was looking in TGCA how the survival is affected ..


text:
Visualising gene expression across cell type and conditions in one plot, in Single Cell Sequencing data
I am new to single cell sequencing data analysis but I have basic programming background in R and python. I want to be able to plot differential expression of two genes, Gene1 and Gene2, across three cell types and across three conditions on one plot. I have followed this tutorial to perform QC and various differential analyses using Seurat package on my data up till now.
My data set is one object, containing information about the three cell types under three various conditions. The metadata slot of my data set contains information about my cell types as well as the conditions under which they are tested. Using the following DotPlot commands I am able to generate separate plots of gene expression with respect to cell type and with respect to condition:
Seurat::DotPlot(object = object, genes.plot = c("Gene1", "Gene2"), plot.legend = TRUE, group.by = "cell_types")

Seurat::DotPlot(object = object,
        genes.plot = c("Gene1", "Gene2"), plot.legend = TRUE, group.by = "tissue_types" )

For elegant representation of my data, I want to be able to have the information of the following plots into one plot.
I will be very grateful for any hint on how I can proceed to achieve an elegant presentation as well as information about other possible ways to graphically represent gene expression data for publication or presentation purposes.
Thank you in advance for your help.

Try using SplitDotPlotGG, this will allow you to color the dots by some additional variable (eg condition).
In your case you would want to do something like this:
library(Seurat)

SplitDotPlotGG(object = object,
               genes.plot = c("Gene1", "Gene2"),
               plot.legend = TRUE,
               group.by = "cell_types",
               grouping.var = "tissue_types")

This would give a result similar to below (an example from https://satijalab.org/seurat/immune_alignment.html).

Thank you very much, it works. However the plot does not appear (plotting space on Rstudio remains blank) when I set plot.legend = TRUE, but the plot does appear without a legend  when I set plot.legend = FALSE. How can I proceed to fix this? @TimStuart

You probably need to make the plotting device bigger