text:
How I move these files based on their names in R or terminal

EDITED

I have list of my samples in first column and corresponding file name for my .vcf files in second column of a file name clin_name.txt like below. For example name of one vcf is LP6008031-DNA_A04_vs_LP6008032-DNA_D03.snp.pass.vcf
OC/AH/120/E/T3  LP6005409-DNA_F01
OC/AH/126/E/T1  LP2000325-DNA_A01
OC/AH/136/E/T3  LP6005409-DNA_E02

This is the link of my file https://www.dropbox.com/s/mvzf47feiqe7x69/clin_name.txt.txt?dl=0
In the second file name pre_post.list I have information of pre or post treatment for each sample like this
pre-treated post-treated
OC/AH/120   OC/AH/065
OC/AH/126   OC/AH/069
OC/AH/136   OC/AH/074
OC/AH/160   OC/AH/084

This is link of this https://www.dropbox.com/s/ikrl1h8yoc0sngx/pre_post.list.xlsx?dl=0
As you suggested I typed
cat pre_post.list | awk '{print $2}' | sed 's@.*/@@'  > post_list.txt
cat pre_post.list | awk '{print $1}' | sed 's@.*/@@'  > pre_list.txt

Then I have
pre-treated
120
126
136
160
182

post-treated
065
069
074
084

in I typed ./test.sh in which I wrote
for item in  $(cat  pre_list.txt | grep -Eo '[0-9]{1,4}');
    do if grep -q $item clin_name.txt ;
    then out=$(grep $item clin_name.txt | awk  '{print $1}');
         cp *$out* /temp/hgig/fi1d18/TRG45/snp/
    else echo post_list $item;  fi;  done;

That returns
cp: cannot stat `*OC/AH/069/S/T3*': No such file or directory
[fi1d18@cyan01 snp]$

Or I run
for item in  $(cat  post_list.txt | grep -Eo '[0-9]{1,4}');
    do if grep -q $item clin_name.txt ;
    then out=$(grep $item clin_name.txt | awk  '{print $1}');
         cp *$out* /temp/hgig/fi1d18/TRG45/snp/pre/
    else cp *$out* /temp/hgig/fi1d18/TRG45/snp/post/;  fi;  done;

That returns
cp: cannot stat `OC/ST/026/S/T1*': No such file or directory
[fi1d18@cyan01 snp]$

I really don't know how to solve this :(

This is more of a job for snakemake or similar. You'd end up following the same logic in R if you really prefer that. Show us what you've tried so far, since there's no reason someone should write a custom program for you.

Assuming the names of all the files from one group are in file named list.
cat pre_post.list | awk '{print $2}' | sed 's@.*/@@'  > post_list.txt
cat pre_post.list | awk '{print $1}' | sed 's@.*/@@'  > pre_list.txt

Now you have two number list
    for item in  $(cat  pre_list.txt | grep -Eo '[0-9]{1,4}');
    do if grep -q $item clin_name.txt ;
    then out=$(grep $item clin_name.txt | awk  '{print $1}');
         cp *$out* /dir/for_pre_files/
    else echo post_list $item;  fi;  done

for post
    for item in  $(cat  post_list.txt | grep -Eo '[0-9]{1,4}');
    do if grep -q $item clin_name.txt ;
    then out=$(grep $item clin_name.txt | awk  '{print $1}');
         cp *$out* /dir/for_post_file/
    else echo Prelist $item;  fi;  done

or you could do this if one that isn't in post belongs in pre
    for item in  $(cat  post_list.txt | grep -Eo '[0-9]{1,4}');
    do if grep -q $item clin_name.txt ;
    then out=$(grep $item clin_name.txt | awk  '{print $1}');
         cp *$out* /dir/for_post_file/
    else cp *$out* /dir/for_pre_file/;  fi;  done

Sorry I tried your kindly code but did not work and I got this error in very  early stage  [fi1d18@cyan01 ~]$ cat pre_post.list | awk '{print $2}' | sed 's@.*/@@'  > post_list.txt
cat: pre_post.list: No such file or directory

You have to change the names in my bash code to reflect the names of the actual files on my computer

Thank you I am getting this the for loop  grep: clin_name.txt: No such file or directory
post_list 069 this being repeated for all of samples

You need to replace all the file names in the script with your file names.

Sorry what is "clin_name.txt" here? I mean where clin_name.txt comes from?  If I know what clin_name.txt is, I might be able to adopt that to my case Thank you

file that contains all OC/AH/120/E/T3  LP6005409-DNA_F01
OC/AH/126/E/T1  LP2000325-DNA_A01
OC/AH/136/E/T3  LP6005409-DNA_E02

Hello @Michael Hearn I tried what you kindly wrote for me edited in my main post but I am getting errors like cp: missing destination file operand after `/temp/hgig/fi1d18/TRG45/snp/pre/'


text:
Uploading files from local computer or web in rstudio cloud
I am using RStudio Cloud for running velocyto R package. So I have tried to import velocyto file from http://pklab.med.harvard.edu/velocyto/mouseBM/SCG71.loom
I tried to upload that like so

But saying
> SCG71 <- read.table("/cloud/project/SCG71.loom", header=TRUE, quote="\"")
Error in make.names(col.names, unique = TRUE) :
  invalid multibyte string at '<89>HDF'
In addition: There were 19 warnings (use warnings() to see them)

Can you help me in very initial of velocyto R?

This seems to be an encoding problem (certain characters not being recognized), see this, this and also this. Probably the file "/cloud/project/SCG71.loom" is encoded with a format that is not included in your "locale".
Check the output of Sys.getlocale() and see if your input file's encoding is listed. If not, you will need to try something like Sys.setlocale("LC_CTYPE", "en_US.UTF-8") or Sys.setlocale("LC_ALL", "English") (see the third link for details).

Thank you, says > Sys.getlocale()
[1] "LC_CTYPE=C.UTF-8;LC_NUMERIC=C;LC_TIME=C.UTF-8;LC_COLLATE=C.UTF-8;LC_MONETARY=C.UTF-8;LC_MESSAGES=C.UTF-8;LC_PAPER=C.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=C.UTF-8;LC_IDENTIFICATION=C"

Did you try Sys.setlocale("LC_CTYPE", "en_US.UTF-8")

Thank you, > Sys.setlocale("LC_CTYPE", "en_US.UTF-8")
[1] "en_US.UTF-8"

But are you getting the same error after setting the locale with Sys.setlocale()? Please don't just run commands and report the output but try to make sense of these. There is quite some information and workarounds to problems like yours in the links I provided.

Velocyto.R
I'm not sure that your issues is from RStudio cloud but more about how you are trying to open the file. The file that you are trying to import is not a velocyto file but a .loom file that can be open using the velocyto.R package.
I'm not very familiar with this kind of of .loom file but on the velocyto.R tutorial, it is specified to do the following:
library(velocyto.R)
ldat <- read.loom.matrices("SCG71.loom")

And you have the following file structures:
> str(ldat)
List of 3
 $ spliced  :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. ..@ i       : int [1:2822201] 74 130 178 192 202 223 225 258 353 418 ...
  .. ..@ p       : int [1:6668] 0 464 786 1469 2502 3290 3760 4974 5388 5659 ...
  .. ..@ Dim     : int [1:2] 24421 6667
  .. ..@ Dimnames:List of 2
  .. .. ..$ : chr [1:24421] "Xkr4" "Rp1" "Sox17" "Mrpl15" ...
  .. .. ..$ : chr [1:6667] "SCG71:AAACGATCCTACTTGT" "SCG71:AAACACTACGTAGAAG" "SCG71:AAACACGGCATCGCAG" "SCG71:AAACAAACAGCAGAAC" ...
  .. ..@ x       : num [1:2822201] 1 1 1 1 1 1 1 1 1 4 ...
  .. ..@ factors : list()
 $ unspliced:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. ..@ i       : int [1:528048] 61 100 339 352 377 445 465 528 617 758 ...
  .. ..@ p       : int [1:6668] 0 115 188 283 431 580 696 994 1127 1232 ...
  .. ..@ Dim     : int [1:2] 24421 6667
  .. ..@ Dimnames:List of 2
  .. .. ..$ : chr [1:24421] "Xkr4" "Rp1" "Sox17" "Mrpl15" ...
  .. .. ..$ : chr [1:6667] "SCG71:AAACGATCCTACTTGT" "SCG71:AAACACTACGTAGAAG" "SCG71:AAACACGGCATCGCAG" "SCG71:AAACAAACAGCAGAAC" ...
  .. ..@ x       : num [1:528048] 1 3 1 1 1 1 1 1 1 1 ...
  .. ..@ factors : list()
 $ ambiguous:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. ..@ i       : int [1:47846] 11915 13649 13830 15638 18068 19585 23490 24345 8142 8767 ...
  .. ..@ p       : int [1:6668] 0 8 13 25 42 56 64 87 97 100 ...
  .. ..@ Dim     : int [1:2] 24421 6667
  .. ..@ Dimnames:List of 2
  .. .. ..$ : chr [1:24421] "Xkr4" "Rp1" "Sox17" "Mrpl15" ...
  .. .. ..$ : chr [1:6667] "SCG71:AAACGATCCTACTTGT" "SCG71:AAACACTACGTAGAAG" "SCG71:AAACACGGCATCGCAG" "SCG71:AAACAAACAGCAGAAC" ...
  .. ..@ x       : num [1:47846] 2 1 1 1 1 2 1 2 1 1 ...
  .. ..@ factors : list()

Alternative 1 - LoomExperiment
On bioconductor, you can also find LoomExperiment. On my computer, once I get it install, I was able to load your SCG71.loom file using:
library(LoomExperiment)
scle <- import("SCG71.loom",type = "LoomExperiment", rownames_attr = NULL, colnames_attr = NULL)

Then, you can find the output here:
> head(scle)
class: LoomExperiment
dim: 6 6667
metadata(4): chunks creation_date velocyto.__version__ velocyto.logic
assays(4): matrix ambiguous spliced unspliced
rownames: NULL
rowData names(6): Accession Chromosome ... Start Strand
colnames: NULL
colData names(1): CellID
rowGraphs(0): NULL
colGraphs(0): NULL

You can then load some data
> scle$CellID[c(1:5)]
[1] "SCG71:AAACGATCCTACTTGT" "SCG71:AAACACTACGTAGAAG" "SCG71:AAACACGGCATCGCAG"
[4] "SCG71:AAACAAACAGCAGAAC" "SCG71:AAACAAACTTTGGGAG"

After that, I don't know what to do with that but I think your issue of importing a datafile is because you don't have the right package to open it.
Alternative - using LoomR
Alternatively, you can also install the package loomR and then enter:
library(loomR)
lfile = connect(filename = "SCG71.loom", mode = "r+")

And you get the output:
> lfile
Class: loom
Filename: /***my_directory***/SCG71.loom
Access type: H5F_ACC_RDWR
Attributes: creation_date, chunks, velocyto.__version__, velocyto.logic, version
Listing:
       name    obj_type dataset.dims dataset.type_class
  col_attrs   H5I_GROUP         <NA>               <NA>
 col_graphs   H5I_GROUP         <NA>               <NA>
     layers   H5I_GROUP         <NA>               <NA>
     matrix H5I_DATASET 6667 x 24421          H5T_FLOAT
  row_attrs   H5I_GROUP         <NA>               <NA>
 row_graphs   H5I_GROUP         <NA>               <NA>

Then, for manipulating these files, you should refer to their tutorial (it's beyond my skills and the scope of your question).
Hope it helps you to figure the solution of your question

source: Stack Exchange

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

sort: votes


text:
Putting cells in Seurat dimplot selectively
I am going to adjust Seurat dimplot in a way avoiding some cells so both my dimplot and heatmap look nice
If you please consider this picture, you would see some cells are far from the clusters so I want to avoid them in dimplot and of course for heatmap (coming from finding markers).
[![enter image description here][1]][1]
Can you suggest me a piece of code doing that?
Thanks

Answered on GitHub here: https://github.com/satijalab/seurat/issues/2338

If you're using a GUI you could select the cells interactively:
plot <- DimPlot(seurat_obj, reduction = "umap")

Then select the cells by clicking around them
select.cells <- CellSelector(plot = plot)
Idents(seurat_obj, cells = select.cells) <- "SubCells"

and subset based on these cells
sub_cells <- WhichCells(seurat_obj, idents = "SubCells")
sub_obj <- subset(seurat_obj, cells = sub_cells)

You could then use subobject in your downstream dimplot.
See here for some of the interactive options. I guess you could outright subset on UMAP coordinates, but that may not always be easy if your boundary isn't linear.

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


text:
Intersecting a data frame based on the multiple columns
I have 2 data frames
> head(df1)
   Chromosome Start_Position End_Position       Gene
1:      chr12       46254706     46254706       ARID2
2:      chr17        7579849      7579849        TP53
3:       chr9       21974777     21974780      CDKN2A
4:       chr1       27105930     27105930      ARID1A
5:       chr1      201982973    201982973        ELF3
6:       chr3       89480459     89480459       EPHA3
> head(df2)
  Chromosome Start_Position End_Position
1          1          10583       863511
2          1       12841835     12854479
3          1       17231299     17232877
4          1      120380739    142540462
5          1      142540484    142684781
6          1      142685823    143462598
>

How I can intersect genes from df1 which are in the range of df2? Because both data frames have chromosome number, start and end but I don't know how to identify genes from df1 based on df2
I have tried
df1[df1$Start_Position%in%df2$Start_Position,]

But I obtain one gene
I also tried
> new_dataset <- df1 %>% right_join(df2, by=c("Chromosome","Start_Position"))
Error: Can't join on 'Start_Position' x 'Start_Position' because of incompatible types (character / integer)
>

reformat df2 so that you have chr1 is instead of 1, then you basically have two bed files.  Then use bedtools intersect.  https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html

In R:
> write.table(df1, 'df1.txt', quote=F, sep="\t", row.names=F, col.names=F)
> write.table(df2, 'df2.txt', quote=F, sep="\t", row.names=F, col.names=F)

Convert both files to sorted BED:
$ sort-bed df1.txt > df1.bed
$ awk -v FS="\t" -v OFS="\t" '{ print "chr"$1 }' df2.txt | sort-bed - > df2.bed

Then map df1 gene names to df2 coordinates:
$ bedmap --echo --echo-map-id-uniq --delim '\t' df2.bed df1.bed > answer.bed

To bring back into R, e.g.:
> dfMap <- read.table('answer.bed', sep='\t', col.names=c('Chromosome', 'Start_Position', 'End_Position', 'Genes'), colClasses=c("factor", "numeric", "numeric", "character"))


text:
Summing up one column entries for each sample
I have a dataframe with 27 samples and for each of them I have different number of values
I want to sum up the values for each sample, leaving each sample with a single value
> dput(d[1:10,c(1,8)])
structure(list(SampleID = c("LP6008460-DNA_A04", "LP6008460-DNA_A04",
"s30", "s30", "s67", "s15", "s67", "s67", "s67", "s30"), X.1 = c(8.17407e-05,
1.40067e-06, 3.86017e-05, 9.50247e-05, 7.08333e-07, 3.19133e-05,
3.1602e-05, 5.00513e-05, 3.29613e-05, 1.003e-06)), row.names = c(NA,
10L), class = "data.frame")
>

I started to do that manually but that was a pain and inaccurate

I have answered this question, but it seems to be more suitable for stack overflow since it's just a programming question without any bioinformatics in. And it's also best to at least pretend you've made an effort to answer the question before, since this is a very common question and already has many answers.

It's worth learning how to use dplyr for this sort of thing:
dat %>%
    dplyr::group_by(SampleID) %>%
    dplyr::summarise(sample_mean = mean(X.1))

# A tibble: 4 x 2
  SampleID          sample_mean
  <chr>                   <dbl>
1 LP6008460-DNA_A04   0.0000416
2 s15                 0.0000319
3 s30                 0.0000449
4 s67                 0.0000288

Sorry for each sample should be sum not mean

I think you can make that tiny change yourself, @Exhausted

Thank you @RamRS I replaced mean with sum myself but did not make sense

What do you mean by "did not make sense"?


text:
move files of a different group to separate directory

EDITED

I have list of my samples in first column and corresponding file name for my .vcf files in second column of a file name clin_name.txt like below. For example name of one vcf is LP6008031-DNA_A04_vs_LP6008032-DNA_D03.snp.pass.vcf
OC/AH/120/E/T3  LP6005409-DNA_F01
OC/AH/126/E/T1  LP2000325-DNA_A01
OC/AH/136/E/T3  LP6005409-DNA_E02

This is the link of my file https://www.dropbox.com/s/mvzf47feiqe7x69/clin_name.txt.txt?dl=0
In the second file name pre_post.list I have information of pre or post treatment for each sample like this
pre-treated post-treated
OC/AH/120   OC/AH/065
OC/AH/126   OC/AH/069
OC/AH/136   OC/AH/074
OC/AH/160   OC/AH/084

This is link of this https://www.dropbox.com/s/nhv3odjkdjdp7hz/pre_post.list.txt?dl=0
For moving different .vcf files based on the pre or post treatment to separate folders I typed
cat pre_post.list | awk '{print $2}' | sed 's@.*/@@'  > post_list.txt
cat pre_post.list | awk '{print $1}' | sed 's@.*/@@'  > pre_list.txt

Then I have
pre-treated
120
126
136
160
182

post-treated
065
069
074
084

I then typed ./test.sh in which I wrote
for item in  $(cat  pre_list.txt | grep -Eo '[0-9]{1,4}');
    do if grep -q $item clin_name.txt ;
    then out=$(grep $item clin_name.txt | awk  '{print $1}');
         cp *$out* /temp/hgig/fi1d18/TRG45/snp/
    else echo post_list $item;  fi;  done;

That returns
cp: cannot stat `*OC/AH/069/S/T3*': No such file or directory
[fi1d18@cyan01 snp]$

Or I run
for item in  $(cat  post_list.txt | grep -Eo '[0-9]{1,4}');
    do if grep -q $item clin_name.txt ;
    then out=$(grep $item clin_name.txt | awk  '{print $1}');
         cp *$out* /temp/hgig/fi1d18/TRG45/snp/pre/
    else cp *$out* /temp/hgig/fi1d18/TRG45/snp/post/;  fi;  done;

That returns
cp: cannot stat `OC/ST/026/S/T1*': No such file or directory
[fi1d18@cyan01 snp]$

I really don't know how to solve this :(

You have quite a few bad shell practices here (never do for i in $(cat file), use while read i; do ... ;done < file instead; you don't need cat file | grep, you can just do grep file), but those shouldn't be causing issues here. Why are you expecting a file to match OC/ST/026/S/T1? Isn't your file name the other column?

Also, are you really cating an excel file? Don't you have it saved as text?

And what is your sample name? For example, you have OC/AH/120/E/T3 in clin_name.txt but OCCAMS/AH/120 in pre_post.list.xlsx. The sample names don't match. All of the names in the pre_post.list.xlsx file start with OCCAMS but OCCAMS doesn't exist in any of the other files.

Thank you @terdon for comments. I changed excel to txt also OCCAMS changed to OC. For example OC/AH/120 is patient key saying pre treated which matches with LP6005409-DNA_F01 in .vcf file. potentially i want to move vcf files based on pre or post treated to separate folders

You need to bring together the information in the two tables, I prefer R for such actions. The common columns between the two tables do not quite contain the same information hence additional data cleaning was required. Moreover the table with treatment information had to be changed to long format with melt. And I have filtered the merged_data so as to conation only "pretreatment" info but you can add a if/else statement to the for loop to control how/where you will move your files.
library(data.table)

clin_name <- fread("clin_name.txt", header = TRUE)
clin_name$sample_id <- sub("/E/T.", "", clin_name$sample_id)

pre_post <- fread("pre_post.list", header = TRUE)
print(pre_post)
pre_post <- melt(pre_post,
     measure.vars = c("pre-treated", "post-treated"),
     variable.name = "treatment",
     value.name = "sample_id")

> print(pre_post)
      treatment sample_id
1:  pre-treated OC/AH/120
2:  pre-treated OC/AH/126
3:  pre-treated OC/AH/136
4:  pre-treated OC/AH/160
5: post-treated OC/AH/065
6: post-treated OC/AH/069
7: post-treated OC/AH/074
8: post-treated OC/AH/084

merged_data <- merge(clin_name, pre_post)
print(merged_data)

pre <- merged_data[merged_data$treatment == "pre-treated", "file_name"]

> print(merged_data)
   sample_id         file_name   treatment
1: OC/AH/120 LP6005409-DNA_F01 pre-treated
2: OC/AH/126 LP2000325-DNA_A01 pre-treated
3: OC/AH/136 LP6005409-DNA_E02 pre-treated

# invoking a terminal command to move each file to a directory called "pre"
for (i in pre$file_name) {
  command <- paste0("mv ", i, " pre/")
  system(command)
}

Thank merged_data returns nothing with zero dimension

merge expects shared column names in the two tables to be merged, if this is not the case you will need to explicitly state the columns to be used when merging, i.e. by using merge(table1, table2, by.x = "column_name_in_table1", by.y = "column_name_in_table2). Line 4 in my code is creating a column called sample_id in clin_name. This column is also present in pre_post (after "melting").

sample_id in clin_name is all NA

You should find a way to create a sample_id column in clin_table and this column should have the same sample encoding as pre_post. The last 5 or so characters of your sample ids in the first table should be removed to make them same. Be careful with the column names, your first table has none.


text:
Attempting to generate a heat map like in this figure
I have a data frame
Response TO CHEMOTHERAPY    Gene    CCF
Responders                  POLQ    1
Responders                 CDKN2A   1
Responders                  TP53    1
Responders              CCDC102B    0.906729075
Non-responders           CCDC102B   1
Non-responders           CCDC102B   1
Non-responders           CCDC102B   1
Non-responders             TSHZ3    1
Non-responders             PIK3CA   1
Responders                 SLIT2    1
Non-responders              SLIT2   1

I have two groups of patients; Responders to chemotherapy and non-responders. For each group I have calculated cancer cell faction (CCF) for a fixed set of genes so that an individual gene has several reads in each group
If CCF = 1 the gene clonal and if CCF < 1 gene is sub clonal, for a given gene for each group, I therefore would have the number of clonal and sub clonal for this gene, for example 2 clonal 7 sub clonal which would be 2/7 something like the below heat map

By your kindly code I produced this

Which is not well informative
This is another potential example heat map

Which show the clonal frequency for each group
By the way thanks for your kindness
I don't know how to use dput() but this is my complete data
https://www.dropbox.com/s/ol7m4nlkuyuy85a/c1.txt?dl=0

What have you tried? You could use geom_rect after setting x and y aes to Sample and Gene.

There is multiple examples of how to do heatmap with ggplot2 on SO, for example, you can start by taking a look here: https://stackoverflow.com/questions/8406394/how-to-produce-a-heatmap-with-ggplot2/8406774#8406774 or https://stackoverflow.com/questions/14290364/heatmap-with-values-ggplot2 or https://stackoverflow.com/questions/60162198/customizing-heatmap-appearance-using-ggplot-in-r/60163171#60163171 ... Try to reproduce those examples and then you will know how to deal with your data

Given the lack of effort from the OP, probably I shouldn't have wasted my time on this question with a lot of problems:
i) The question is not reproducible, the OP should have supplied their data with dput().
ii) The data as pasted above does not quite make sense to me, why would you have three readings from the same gene and even in such a case how should these replicates be represented in a heatmap?
iii) The question is not specific enough and should not go like "how can I generate a heatmap" and should be more like "in an attempt to generate a heatmap like in this figure I have tried this piece of code but ..."
Having said these, and after downvoting the question, here is a solution leveraging the ComplexHeatmap package:
library(data.table)
library(ComplexHeatmap)

my_data <-
"
Sample  Gene    CCF
Sample1 POLQ    1
Sample1 CDKN2A  1
Sample1 TP53    1
Sample1 CCDC102B    0.906729075
Sample1 CCDC102B    1
Sample1 CCDC102B    1
Sample1 CCDC102B    1
Sample2 TSHZ3   1
Sample2 PIK3CA  1
Sample2 SLIT2   1
Sample2 SLIT2   1
"

my_data <- fread(my_data)

my_data <- dcast.data.table(my_data,
                            formula = Gene ~ Sample,
                            value.var = "CCF",
                            fun=mean)

# dcast generates NA values as the provided data does not have values for each gene
my_data[is.na(my_data)] <- 0

# Heatmap() expects a matrix input
my_matrix <- as.matrix(my_data[,c("Sample1", "Sample2")])
rownames(my_matrix) <- my_data$Gene

Heatmap(my_matrix,
        cell_fun = function(j, i, x, y, width, height, fill) {
          grid.text(my_matrix[i, j], x, y, gp = gpar(fontsize = 10))})

This is how the data look after dcast, in which multiple entries per gene (my point ii above) is summarized with mean():
> head(my_matrix)
           Sample1 Sample2
CCDC102B 0.9766823       0
CDKN2A   1.0000000       0
PIK3CA   0.0000000       1
POLQ     1.0000000       0
SLIT2    0.0000000       1
TP53     1.0000000       0

Here is the resulting heatmap:

Sorry please don't close my question, I have edited that and I am certain I know what I want and what is the final goal


text:
Error in ruuning a package in R
I want to install Signature Tools Lib R package
I typed
> devtools::install()
Error: Could not find package root, is your working directory inside a package?
>

Packages does need this list of sotware
R package dependencies:
    VariantAnnotation,
    BSgenome.Hsapiens.UCSC.hg38,
    BSgenome.Hsapiens.1000genomes.hs37d5,
    BSgenome.Mmusculus.UCSC.mm10,
    BSgenome.Cfamiliaris.UCSC.canFam3,
    SummarizedExperiment,
    BiocGenerics,
    GenomeInfoDb,
    NMF,
    foreach,
    doParallel,
    doMC,
    lpSolve,
    ggplot2,
    methods,
    cluster,
    stats,
    NNLM,
    nnls,
    GenSA,
    gmp,
    plyr,
    RCircos,
    scales,
    GenomicRanges,
    IRanges,
    BSgenome

I have installed all of these one by one but again when trying things like the below, I get error
Error in tabToSNVcatalogue(subs = tmpSNVtab, genome.v = "hg19") :
  could not find function "tabToSNVcatalogue"
> devtools::test()
Error: Could not find package root, is your working directory inside a package?
>
> library(signature.tools.lib)
Error in library(signature.tools.lib) :
  there is no package called ‘signature.tools.lib’
>

First, we would need to know which is the package you want to install. Second this is related to Bioinformatics? it doesn't seem so, could you edit the question to include this details?

If you do ?devtools::install, you will see that the install() function from the devtools package, not surprisingly, requires a package to be specified:
install(pkg = ".", reload = TRUE, quick = FALSE, build = !quick,
  args = getOption("devtools.install.args"), quiet = FALSE,
  dependencies = NA, upgrade = "ask", build_vignettes = FALSE,
  keep_source = getOption("keep.source.pkgs"), force = FALSE, ...)

So it is absolutely normal that devtools::install() does not install anything. Just specify the packages(s) with the pkg argument.
And this is how you install Signature Tools:
install.packages("remotes")
remotes::install_github("Nik-Zainal-Group/signature.tools.lib")

Thank you but your solution gives ERROR: dependency ‘NNLM’ is not available for package ‘signature.tools.lib’

removing ‘/Library/Frameworks/R.framework/Versions/3.6/Resources/library/signature.tools.lib’
Error: Failed to install 'signature.tools.lib' from GitHub:
(converted from warning) installation of package ‘/var/folders/wk/qwgr9yqs6ts0czx__y_w8gch0000gq/T//RtmpA8PudK/filec5bcb188c2b/signature.tools.lib_0.0.0.9000.tar.gz’ had non-zero exit status


text:
Getting genes specially up or down regulated
I have 6 RNA-seq samples like this
4 patients (005, 036, 121, 013)

I have 3 tumour samples and 3 cancer models (organoid)

This is PCA of log transformed data by DESeq2 (tumour Vs. organoid)

In this paper https://www.nature.com/articles/s41467-018-05190-9
They are saying

Finally, for each patient, the genes that are specifically up- or
  down regulated in each patient’s organoid (1 vs. all comparison, p ≤
  0.05 and abs (log2-FoldChange) ≤ 1) and that are among the top-50 genes with highest base mean expression were selected.

They have tumour and organoid like my case; My confusion is what they mean by "each patient’s organoid (1 vs. all comparison"? Because like them I need up or down regulated genes in each patient's organoid but I can not figure out I should compare each organoid versus what?
The ultimate goal would finding a signature of genes specific to each organoid as each organoid originates from tumour of a patient by these signatures one would be able to classify cancer to sub types

Nice results! It will be very easy to find the responsible genes with the tightness of those clusters You perform eg. a linear regression (GLM or machine learning) and look at the regression weights (genes). The results are so clean any method of classification will identify them. Just to mention I don't do cancer (never have, never will), but those results look cool.

Thanks a lot for your positive and encouraging comments. The assumption is that the expression pattern in each organoid should be very close and similar to the tumour from which this organoid originates. More close better an organoid would be a model of that tumour. But in my PCA we are seeing tumours are in one way and organoids are another way. My expectation was for example tumour 005 and organoid 005 placed close to each other. However I will need a signature of each oganoid as a model of cancer. You meant I take genes as independent variables and a specific organoid as response variable?

@MichaelG. Could you please elaborate on how you would use GLM to do this? Or could you send me the names of some of the packages you have used to do this? I'd like to try it out myself.. Thanks.

@FereshTeh, just came across this rather old post. Regarding "My expectation was for example tumour 005 and organoid 005 placed close to each other" part of your comment above, please see the last paragraph of the "RNA-seq analysis" section of the methods from the paper that you shared. They kick out genes that they think arise from the "artificial environment" (their own words) organoids are prone to. You would get a better looking PCA plot if you would follow their practice, however, organoid biology would still be very very (the nuance is arguable) different from tumor biology.


text:
Variant calling without matched normal sample
I have WGS .bam files for 3 patients (tumour and its matched derived model namely organoid) but I don't matched normal sample. If I call variants of each patients (tumour and its matched organoid), how I can use read counts at germline heterozygous positions estimated by GATK 3.2-2 to compensate for the absence of matched normal sample? I heard people use dbsnp VCF instead of the matched normal small variant VCF. But, I don't know start from where? which GATK function does that?

I should mention, Calling SNV and indel in many tools returns vcf but
  I called copy number by varscan that did not return .vcf output so I
  am not sure what to do for CNV

Any helps please?

If you just want to filter out calls present in dbSNP then use:
java -jar GenomeAnalysisTK.jar \
   -T SelectVariants \
   -R reference.fasta \
   -V patient.vcf \
   --discordance dbSNP.vcf \
   -o patient.filtered.vcf

--discordance will produce calls not present in dbSNP.

Thanks a lot @Devon, Calling SNV and indel in many tools returns vcf but I called copy number by varscan that did not return .vcf output, is this command is applicable yet?

It's usually best to mention things like that in your post... I presume this will still work, but I can't say I've ever tried.

https://software.broadinstitute.org/gatk/documentation/article?id=11682

I heard GenomeAnalysisTK.jar has been replaced. Do you know by which function we can filter a vcf file based on dbSNP?

Here is a clue by scatngs
https://github.com/cancerit/ascatNgs/wiki/Human-reference-files-from-1000-genomes-VCFs
and already calculated file is here
ftp://ftp.sanger.ac.uk/pub/cancer/support-files/CPIB/ascatNgs/Human/GRCh37/
Although installing scatngs itself is a big burden


text:
Different contrasts in DESeq2
I have a treatment and control in two time points like this
 > design
                                 X condition time
    1       CTRL_24_hrs_replicate1   control   24
    2       CTRL_24_hrs_replicate2   control   24
    3       CTRL_24_hrs_replicate3   control   24
    4  treatment_24_hrs_replicate1        t   24
    5  treatment_24_hrs_replicate2        t   24
    6  treatment_24_hrs_replicate3        t   24
    7       CTRL_48_hrs_replicate1   control   48
    8       CTRL_48_hrs_replicate2   control   48
    9       CTRL_48_hrs_replicate3   control   48
    10 treatment_48_hrs_replicate1        t   48
    11 treatment_48_hrs_replicate2        t   48
    12 treatment_48_hrs_replicate3        t   48
    >

I want to test between between treatment and control considering time point 2 hours to 4 hours
I have done like this
dds <- DESeqDataSetFromMatrix(countData=a,colData=design, design=~time + condition + time:condition)

ddsTC <- DESeq(dds, test="LRT", reduced = ~ time + time:condition)

> resultsNames(ddsTC)
[1] "Intercept"                "condition_t_vs_control"
[3] "time"                     "conditiont.time"

What is the difference of the results of these codes
results(ddsTC, name="condition_t_vs_control", test="Wald")

versus
results(ddsTC, name="conditiont.time", test="Wald")

Actually I want to know what condition_t_vs_control gives and what conditiont.time gives
Thank you so much for any intuition

condition_t_vs_control gives the effect of condition, conditiont.time is the interaction of condition and time. If you wanted to test the effect of time, use name="time".

Thank you, using which contrast I can get the difference of IT versus control but considering time goes from 24 hours to 48 hours?

Just condition_IT_vs_control, since time is compensated for by fitting all coefficients at the same time (let me know if that phrasing is confusing).

Sorry @Devon Ryan I get  error for ddsTC <- DESeq(dds, test="LRT", reduced = ~ time + time:condition)                                                                                error in nbinomLRT(object, full = full, reduced = reduced, quiet = quiet,  :
less than one degree of freedom, perhaps full and reduced models are not in the correct order  what should I do please ?

Make sure to specify the full model, I don't know what that defaults to.

Thank you so much for your answer, I think this is my full model                                              dds <- DESeqDataSetFromMatrix(countData=a,colData=design, design=~time + condition + time:condition)   am I right? design=~time + condition + time:condition  When I try to reduce that by removing reduced = ~ time + time:condition I get error

Then convert time to a factor. Perhaps that will solve the problem. You have enough samples that you shouldn’t run into this issue. BTW you can use the default Wald test instead and avoid the error.

Thank you so much says that > ddsTC <- DESeq(dds,reduced = ~ time + time:condition)
Error in DESeq(dds, reduced = ~time + time:condition) :
'reduced' ignored when test='Wald'

Correct, there is no reduced formula with a Wald test, you just extract the appropriate coefficients with lfcShrink()


text:
DEseq2 contrast in two time points
I need help with some clarifications please
I have control versus treatment in two time points like
> y
                             X condition time
1       CTRL_24_hrs_replicate1   control   24
2       CTRL_24_hrs_replicate2   control   24
3       CTRL_24_hrs_replicate3   control   24
4  treatment_24_hrs_replicate1        t   24
5  treatment_24_hrs_replicate2        t   24
6  treatment_24_hrs_replicate3        t   24
7       CTRL_48_hrs_replicate1   control   48
8       CTRL_48_hrs_replicate2   control   48
9       CTRL_48_hrs_replicate3   control   48
10 treatment_48_hrs_replicate1        t   48
11 treatment_48_hrs_replicate2        t   48
12 treatment_48_hrs_replicate3        t   48
>

I want to see the difference of treatment versus control from 24 hours to 48 hours, so this is not just the difference of treatment versus control or difference between 24 versus 48 hours rather the treatment versus control from 24 hours to 48 hours
This is what I have done so far
>  dds <- DESeqDataSetFromMatrix(countData=count,
+                                 colData=design,
+                                 design=~condition, tidy = FALSE)
Warning message:
In DESeqDataSet(se, design = design, ignoreRank) :
  some variables in design formula are characters, converting to factors
>   dds$group <- factor(paste0(dds$time, dds$condition))
>   design(dds) <- ~ group
>   dds <- DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
>   resultsNames(dds)
"Intercept"                    "group_24IT_vs_24control"         "group_48control_vs_24control" "group_48IT_vs_24control"
>

How to find the difference at time point 24 vs 48 in condition ?
I searched and I found a similar thing but I don't know how to adopt that to my data
results(dds, contrast=list(c("time_4_vs_0", "conditionLD.time4")), test="Wald")

I tried this and I got error
> results(dds, contrast=list(c("time_48_vs_24", "conditionIT.time24")), test="Wald")
Error in checkContrast(contrast, resNames) :
  all elements of the contrast as a list of length 2 should be elements of 'resultsNames(object)'
>

Can you help me in adopting this to my data please?

You can't figure out how to adapt this to your data?
http://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#contrasts
Ignore resultsNames.  You don't need them for simple comparisons like you want to do.

Sorry even by ignoring resultsNames, I don't know what to put into this formula. For instance what I should put instead of time_4_vs_0 ?   This gave me error  > results(dds, contrast=list(c("time_48_vs_24", "conditionIT.time24")), test="Wald")
Error in checkContrast(contrast, resNames) :
all elements of the contrast as a list of length 2 should be elements of 'resultsNames(object)'

Where have you looked?  Did you read the link I posted?  Did you look at the entry for deseq2 results?  https://rdrr.io/bioc/DESeq2/man/results.html   Can you not find a single tutorial anywhere which shows you how to use specify the contrast?

you can do:
data$condition_new <- as.factor(paste0(data$condition, " - ", data$time))
dds <- DESeqDataSetFromMatrix(countData=count,
                             colData=design,
                             design=~condition_new
                             , tidy = FALSE)
enter code here
results(dds, contrast=c("condition_new", "control - 24", "condition - 4"))

I hope it helps.

You can find here: http://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#contrasts
the tip as pointed from @swbarnes2


text:
Finding mutations in glycosylation sites
We are going to look at the likely N- and O- glycosylation sites within MUC16 (Q8WXI7 · MUC16_HUMAN)§ by in house long-read DNA sequencing data (PacBio).
Counting the number of tandem repeats is a task which could be done with bioinformatic tools like ad-codeVNTRas shown in the below figure. Another task is identifying mutations within these repeats. My question is: after predicting glycosylation sites using tools like PA2DE within these tandem repeats in MUC16 gene, how could I identify any point mutations in these glycosylation sites?
I googled a lot to get some idea; Shall we do germline analysis in MUC16 glycosylation predicted sites? Is this what we should use to look at glycosylation sites that are ‘rejected’ to understand the sites that can’t be glycosylated (mutated)?
Does anyone have such an experience how to identify mutation in glycosylation sites in MUC16 which known to have variable number tandem repeat (VNTR)?

If you write the exact rules, it's not difficult to write a regex function to find these. How many repeats? What does tandem mean? Is it without spaces, VNTRVNTRVNTR?

Thank you, counting the number of repeat is another task which could be done with bioinformatic tools like ad-codeVNTR. Another task is identifying mutations within these repeats. My question is after predicting glycosylation sites, how could I identify any point mutations in these glycosylation sites.

Please [edit] your question and clarify what you need and what you are working with:  What are your input data? Protein? DNA? Have you performed your own sequencing? Are you looking through published sequences? How do you define mutation, what is your reference? Doesn't code-adVNTR do exactly what you want? The image you show suggests it does, so isn't it already what you need?

Thank you I tried to edit my question. That would DNA long read sequencing done on human cancer cell lines


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:
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:
Fitting a principal curve into a diffusion map
I returned a Seurat object with a diffusion map by
 seurat_object <- RunDiffusion(seurat_object,genes.use = seurat_object@var.genes)

as below

Now, how I could know the order of a cell's projection onto this principal curve representing its predicted progression through the development by fitting a principal curve on my diffusion map?
I know Seurat slot for DM1 and DM2 is inside the seurat@dr
by
fit=principal_curve(as.matrix(seurat_object@dr$dm@cell.embeddings))
plot(fit)
line(fit)
points(fit)

I draw this graph of my cells

But how I know the first and last cells on curve? I mean how I know cells from each cluster is being placed where on the curve? By fit$ord I will obtain orders of cells on the line but how I know the first and end of the line?????

If fit$ord gives you the order of cells on the curve, you should be able to obtain the first and last cells by doing the following:
first_cell_index <- which.min(fit$ord)
last_cell_index <- which.max(fit$ord)

However, I'm not sure this is so informative - if you're looking to get representative gene expression for the start and end of the curve, you should really be looking at the cluster of cells towards the start and end of the curve - for example, by comparing cells in clusters 1 and 3 that you showed in your first plot. A single cell will suffer from dropout and not tell you a lot about the cluster.

Thank you, fit$lambda gives me pseudotime amount; Now I want to have a violin plot for lambda for each cell and each cluster


text:
Variability in genes across platforms
I have a matrix of gene expression for single cell (cells in column and genes in row) in 9 time points with about to 200 cells for each time point. I have also a matrix of gene expression in bulk RNA-seq with the same time point and 2 replications.
I have pooled reads from cells in each time point to have a simulated bulk RNA-seq. now I have two matrices; Real bulk RNA-seq in 9 time points and two replication and another matrix a simulated bulk RNA-seq with 9 time points and two replications. I have plotted a heat map of them to show the pattern of gene expression is similar in both matrices so single cell RNA-seq is good enough. Actually I have used the real bulk RNA-seq in a similar study as a gold standard for checking my single cell RNA-seq. My problem is this the pattern of gene expression in 0 time point is very different in two datasets. I should note that 0 time point comes from Fluidigm and another 8 time points come from icell8. I should be able TO adjust this heat map to look similar enough.
So I thought to show similarity by plotting cumulative fraction of genes detectable in varying fractions of cells (lets say 20%, 40%, ...), segregated by bulk RNA-seq expression quartile. Something like a plot that x axis is percent of cells and y axis is percent of genes but I think I should solve the fundamental problem in a heat map.
I have single cell (matrix with T) and bulk RNA-seq (matrix with R). Whatever I am plotting, time point 0 in single cell RNA-seq seems more yellow and dose not allow a good similarity between single cell and bulk RNA-seq. How I can make zero time point less yellow so another time points become more yellow?

How do you rank each cell by number of genes expressed?

Thank you, actually I am not sure, if I have FPKM values for bulk RNA-seq, I don't   know how to show how many cells express a quartile of genes in bulk RNA-seq

Sorry, I think I misunderstood the question: what has to do the bulk RNA-seq with identifying the variability of genes in individual cells? What do you try to show ? You expect that this graph (however it is made) will help you to find/learn that ...?

You right the question is very ambiguous,  I want to check similarity of two sets of samples

First make sure what you want help with in this question. Then please [edit] the question with more details. Do you want to compare the bulk RNA to the scRNA? Or you want to check the similarity at certain timepoints or between replications?

I'm kind of surprised that your method of comparing data from two totally different sample prep methods worked at all, not that it differs significantly at only one timepoint.

What do you want? A method to normalize between technologies so that your time point 0 looks more similar between scRNA and bulk or plot the similarity  between bulk and scRNA? Also I assume that this data wasn't corrected by batch/technology effect or it is?

I have not done batch correction. I only have pooled 200 cells in each time point to simulate a bulk RNA-seq and I hoped when comparing that with real RNA-seq from the same experiment the data look similar. But likely 0 time point is very different coming from two sequencing methods. after pooling my cells I normalised them with deseq2 like I have done for real bulk RNA-seq. for 8 time points looks good but for zero as this platform gives more read counts. how I manipulate my matrices to make a more similar heat maps or how could I remove difference between 0 time point and another time points

Sorry you are right about different sample prep but when I am calculation Pearson correlation between each time point I am seeing a high correlation between data. If this high correlation is by chance? correlation after log transforming both of datasets


text:
Ordering the cells by the expression of specific genes
I have 209 cells. I want to show the expression of two genes by ordering my cells. Imagine something like a plot on which y axis shows the read counts range from 0 to 10000 and x axis shows the number of cells (there is a link of an example picture that I wish to produce the same). So that two colours bars each colour for a gene.
I have used this code but result is not as I want at all
# Library
 library(tidyverse)

 # Create data
 value1=abs(rnorm(26))*2
 data=data.frame(x=LETTERS[1:26], value1=value1, value2=value1+1+rnorm(26, sd=1) )

 # Reorder data using average?
 data = data %>% rowwise() %>% mutate( mymean = mean(c(value1,value2) )) %>% arrange(mymean) %>% mutate(x=factor(x, x))

 # plot
 ggplot(data) +
   geom_segment( aes(x=x, xend=x, y=value1, yend=value2), color="grey") +
   geom_point( aes(x=x, y=value1), color=rgb(0.2,0.7,0.1,0.5), size=3 ) +
   geom_point( aes(x=x, y=value2), color=rgb(0.7,0.2,0.1,0.5), size=3 ) +
   coord_flip()

This is link of my data and my result picture

And to some data

Could you add an example of what a lollipop is? (a link or something like that?) And what have you tried to make that plot (To avoid suggesting to do the same and help you more specifically)

Thank you, I edited my post. Actually lollipop plot itself dose not matter that can be replaced by any other plots. The problem is that I could not figure out how to plot the same as picture and what I have tried up to now did not answer my question :(

could you further explain how this picture doesn't help you ?. I am still trying to figure out how I can help

Thank you, because on gene seems to be compressed at the end  and does not allow to another gene to be cleared. I don't know how to have another gene on another axis or showing y axis in log10

So, your problem is how to convert an axis to log10? This is answered in Stack Overflow and other venues. In ggplot2 you can't have two axis by design: it is against the grammar of graphics.


text:
How to debug a bash script
Please somebody debug my code
In sub-folders there are 500000 vcf files from which I only need 500 vcf files listed in temp.txt which randomly distributed in sub-folders 10 to 60. The code reads each of 500000 vcf files one by one, if a file was in temp.txt, the code runs
dx run swiss-army-knife -icmd="tabix -h $vcf -R $REGIONS_FILE > $output_file" -iin="$filename_root" -iin="$vcf.*" -iin="/Bulk/DRAGEN\ WGS/Whole\ genome\ variant\ call\ files\ (VCFs)\ (DRAGEN)\ [500k release]/file.bed" --destination "$OUTPUT_FOLDER" --detach --priority high --yes
    done < "$VCF_LIST_FILE"

on that if not, the code goes to next file such a process takes years. I want this code only locates files listed in temp.txt and run the main code on them by iteration.
#!/bin/bash

# Set the working directory to the desired location within DNAnexus

# Number of concurrent processes
N=4

# Iterate from 10 to 60
for folder_number in {10..60}; do
  (
    # Define the input files and parameters for this iteration
    folder_path="/Bulk/DRAGEN\ WGS/Whole\ genome\ variant\ call\ files\ (VCFs)\ (DRAGEN)\ [500k release]/$folder_number"
    VCF_LIST_FILE="temp.txt"
    REGIONS_FILE='file.bed'
    OUTPUT_FOLDER='/out'  # Replace with your desired output folder

    # Set the working directory to the folder for this iteration
    dx cd "$folder_path"

    # Read each VCF file from the list
    while read -r vcf; do
      filename_root=$(basename "$vcf")
      output_file="${filename_root}_extracted.vcf.gz"

      # Use dx run swiss-army-knife to execute tabix on the VCF file
      dx run swiss-army-knife -icmd="tabix -h $vcf -R $REGIONS_FILE > $output_file" -iin="$filename_root" -iin="$vcf.*" -iin="/Bulk/DRAGEN\ WGS/Whole\ genome\ variant\ call\ files\ (VCFs)\ (DRAGEN)\ [500k release]/file.bed" --destination "$OUTPUT_FOLDER" --detach --priority high --yes
    done < "$VCF_LIST_FILE"
  ) &

  # Limit the number of concurrent processes
  ((i=i%N)); ((i++==0)) && wait
done

# Wait for all background jobs to finish
wait

# Finished message
echo "++ Finished extracting genomic coordinates from VCF files"

This is how temp.txt looks. The first two digits of vcf file names defines which sub-folder they are located. For instance, the below vcf files distributed in sub-folders 32, 33 and 34
head temp_new.txt
head: temp_new.txt: No such file or directory
izadia@LV19Y7325V Downloads % head new_temp.txt
11/1195100_24053_0_0.dragen.hard-filtered.vcf.gz
12/1204363_24053_0_0.dragen.hard-filtered.vcf.gz
12/1246629_24053_0_0.dragen.hard-filtered.vcf.gz
13/1311312_24053_0_0.dragen.hard-filtered.vcf.gz
14/1427973_24053_0_0.dragen.hard-filtered.vcf.gz
14/1435043_24053_0_0.dragen.hard-filtered.vcf.gz
15/1558824_24053_0_0.dragen.hard-filtered.vcf.gz
16/1683316_24053_0_0.dragen.hard-filtered.vcf.gz
16/1692707_24053_0_0.dragen.hard-filtered.vcf.gz
17/1749825_24053_0_0.dragen.hard-filtered.vcf.gz
izadia@LV19Y7325V Downloads % tail new_temp.txt
52/5230856_24053_0_0.dragen.hard-filtered.vcf.gz
52/5257892_24053_0_0.dragen.hard-filtered.vcf.gz
54/5453459_24053_0_0.dragen.hard-filtered.vcf.gz
55/5514349_24053_0_0.dragen.hard-filtered.vcf.gz
55/5526010_24053_0_0.dragen.hard-filtered.vcf.gz
57/5718723_24053_0_0.dragen.hard-filtered.vcf.gz
58/5868836_24053_0_0.dragen.hard-filtered.vcf.gz
58/5875882_24053_0_0.dragen.hard-filtered.vcf.gz
58/5889611_24053_0_0.dragen.hard-filtered.vcf.gz
59/5917345_24053_0_0.dragen.hard-filtered.vcf.gz
izadia@LV19Y7325V Downloads %

After running the script I get this error
dxpy.exceptions.DXCLIError: Value provided for input field "in" could not be parsed as file: could not resolve "1195100_24053_0_0.dragen.hard-filtered.vcf.gz" to a name or ID
dxpy.exceptions.DXCLIError: Value provided for input field "in" could not be parsed as file: could not resolve "1204363_24053_0_0.dragen.hard-filtered.vcf.gz" to a name or ID
dxpy.exceptions.DXCLIError: Value provided for input field "in" could not be parsed as file: could not resolve "1246629_24053_0_0.dragen.hard-filtered.vcf.gz" to a name or ID
dxpy.exceptions.DXCLIError: Value provided for input field "in" could not be parsed as file: could not resolve "1311312_24053_0_0.dragen.hard-filtered.vcf.gz" to a name or ID
dxpy.exceptions.DXCLIError: Value provided for input field "in" could not be parsed as file: could not resolve "1427973_24053_0_0.dragen.hard-filtered.vcf.gz" to a name or ID
dxpy.exceptions.DXCLIError: Value provided for input field "in" could not be parsed as file: could not resolve "1435043_24053_0_0.dragen.hard-filtered.vcf.gz" to a name or ID
dxpy.exceptions.DXCLIError: Value provided for input field "in" could not be parsed as file: could not resolve "1558824_24053_0_0.dragen.hard-filtered.vcf.gz" to a name or ID
dxpy.exceptions.DXCLIError: Value provided for input field "in" could not be parsed as file: could not resolve "1683316_24053_0_0.dragen.hard-filtered.vcf.gz" to a name or ID
dxpy.exceptions.DXCLIError: Value provided for input field "in" could not be parsed as file: could not resolve "1692707_24053_0_0.dragen.hard-filtered.vcf.gz" to a name or ID
dxpy.exceptions.DXCLIError: Value provided for input field "in" could not be parsed as file: could not resolve "1749825_24053_0_0.dragen.hard-filtered.vcf.gz" to a name or ID

because the codes starts searching from sub-folder 10 for each of files in new_temp.txt one by one

So what do you want to do? If you don't want to check if the files exist, then don't check if the files exist. Are there any vcf files in your various folders that should not be processed? What is the point of temp,txt? And what does it contain? Full paths? FIle names?

Oh, and your code never checks if a file exists. Just [edit] your question, explain what you really want to do, and maybe we can help.

I edited my main post

And what does temp.txt look like? Please show us an example. We need to know if it is one file name per line, many per line, if it is just the names or the full paths etc.

While you can do it in bash often switching to ie python gives you more debugging options, cleaner code, etc.

Please do not edit your posts as follow up to answers without making the edits extremely clear to follow. Right now, we need to read the answers and your edit history (not question) side by side to understand what's going on and that's a huge pain.

Please don't forget to mark the answer below as "accepted" if this addresses your question (it does look like it).

Thanks for the answer but the problem have not been solved yet so the code still stats checking from sub-folder 10 for the listed files one by one which takes forever

The first two digits of vcf file names defines which sub-folder they are located.

You could use this information to edit temp.txt so it contains the full path:
grep -v ".tbi$" temp.txt | sed 's/^\(.\)\(.\)/\1\2\/\1\2/' > new_temp.txt

I'm assuming here that you don't actually want to run dx run swiss-army-knife on the vcf indexes.
This will turn the example of temp.txt you gave into this:
32/3295094_24053_0_0.dragen.hard-filtered.vcf.gz
33/3309510_24053_0_0.dragen.hard-filtered.vcf.gz
34/3468518_24053_0_0.dragen.hard-filtered.vcf.gz
34/3471960_24053_0_0.dragen.hard-filtered.vcf.gz

You can then simply run a single for loop over that, instead of the nested one you are using now.
eg.
while read vcf
  do
    dx run swiss-army-knife [...]
  done <  new_temp.txt


text:
Target genes for piRNA
Where I can find a database or tool to give me the target genes of PIWI (piRNAs) in human?
I found one but works for worm like piRTarBase
There are a lot for miRNAs around but not for piRNAs

Update
It does indeed appear that while piRBase etc. have human piRNAs, they do not have target information.
Looking at related questions, recommendations are made to map piRNA sequences to the genome to find the targets. One tool recommended is PiPipes, which might assume read input but otherwise might give some ideas.
Original answer
piRBase has this. See "target mRNA" field.

Thank you. Can I use that for human? because all are for non human likely

@Exhausted You are right- how odd. THey have human piRNAs in the db but none of them seem to have targets. I missed that, sorry. I will see if I can figure it out and if


text:
Featureplot by co-expression of some genes
I want to identify MDSCs cells which don't not express HLA-DRs (HLA-DR negative) but express CD14, CD33, and ITGAM
For this, I firstly plotted a feature plot to locate HLA-DR by
cd.features <- list(grep(pattern = '^HLA-DR', x = rownames(x = pbmc_small), value = TRUE))

pbmc_small <- AddModuleScore(object = pbmc_small, features = cd.features, ctrl = 5, name = 'CDFeatures')

FeaturePlot(object = pbmc_small, features = 'CDFeatures1')&
+     scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "RdBu")))

Which returned

For CD14, CD33, and ITGAM I done
> cd.features=list(c("CD14", "CD33",  "ITGAM"))
> pbmc_small <- AddModuleScore(object = pbmc_small, features = cd.features, ctrl = 5, name = 'CDFeatures')
> FeaturePlot(object = pbmc_small, features = 'CDFeatures1')&
+     scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "RdBu")))

>

Which returned

By looking at these two plots I really can not say which part of plot is HLA-DR negative but expressing CD14, CD33, and ITGAM
Do you know a way to do this?
Plotting options don't help because looking at the below plots you can not say which cells are HLA negative and CD14 positive

Will the user who flagged for this post to be closed please add a comment stating why they think this post should be closed? I don't see any problem with it.

My small comment, I'd simply use ggplot2 if I was using R (not hugely helpful). It is a lot more powerful though. @RamRS point: It is a very clear question, nicely presented on a good subject. I can make a good educated guess why someone did this, but its really for meta. Closing questions is a current topic there. It is unlikely the question will be closed - its just a vote. (For the record I upvoted.)

@haci will know this, but will have to get into cellular immunology.

Thinking about it you really want a third axis and a good viewing angle. Plotly might be a better option than featurePlot. Interesting but not my thing.

Very old question, but recently modified, which means that someone might still be looking into this.
I usually prefer to use Ucell as the signatures are calculated within cells, which makes comparisons a bit easier. In addition is allows the formulation of signatures with both positive and negative genes.
hladr_genes <- grep(pattern = '^HLA-DR', x = rownames(x = pbmc_small), value = TRUE)
signatures_list <- list(MDSCs = c("CD14+","CD33+","ITGAM+",paste0(hladr_genes,"-")))

> signatures_list
$MDSCs
[1] "CD14+"     "CD33+"     "ITGAM+"    "HLA-DRA-"  "HLA-DRB5-" "HLA-DRB1-"

# calculate signature
pbmc_small <- AddModuleScore_UCell(pbmc_small,signatures_list,ncores = 8,name = "")
# perform smoothing
pbmc_small <- SmoothKNN(pbmc_small, signature.names = names(signatures_list), reduction="pca")

You need to make a transformed variable to the gene sets of interest. I would rescale the data on a range of 0-1. Then you can make a log2 ratio of the sum(set1)/sum(set2). You could explore other transformations, but this will give contrasting values for high set1 and high set2, but neutral if the expression is similar between the sets.

Your answer could be improved with additional supporting information. Please [edit] to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center.


text:
The correct design for time series experiment
I have a treatment and control in two time points like this
> design
                                  condition time
    1       CTRL_24_hrs_replicate1   control   24
    2       CTRL_24_hrs_replicate2   control   24
    3       CTRL_24_hrs_replicate3   control   24
    4  treatment_24_hrs_replicate1        t   24
    5  treatment_24_hrs_replicate2        t   24
    6  treatment_24_hrs_replicate3        t   24
    7       CTRL_48_hrs_replicate1   control   48
    8       CTRL_48_hrs_replicate2   control   48
    9       CTRL_48_hrs_replicate3   control   48
    10 treatment_48_hrs_replicate1        t   48
    11 treatment_48_hrs_replicate2        t   48
    12 treatment_48_hrs_replicate3        t   48
    >

I want to test between treatment and control considering time point 24 hours to 48 hours
I have done like this
dds <- DESeqDataSetFromMatrix(countData=a,colData=design, design=~time + condition + time:condition)

But at this part I get error, although I am trying different things

> ddsTC <- DESeq(dds, test="LRT", reduced = ~ time + time:condition)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
Error in nbinomLRT(object, full = full, reduced = reduced, quiet = quiet,  :
  less than one degree of freedom, perhaps full and reduced models are not in the correct order

or
> ddsTC <- DESeq(dds, test="LRT", reduced = ~time:condition)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
Error in nbinomLRT(object, full = full, reduced = reduced, quiet = quiet,  :
  less than one degree of freedom, perhaps full and reduced models are not in the correct order
>

I tried these with no error although I am not certain if this makes sense at all while the results of both is the exactly the same
> ddsTC <- DESeq(dds, test="LRT", reduced = ~ time + condition)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
>

> ddsTC <- DESeq(dds, test="LRT", reduced = ~ time)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

Al I want is getting the difference of treatment versus control but considering time goes from 24 hours to 48 hours
Thanks for any help


text:
maftools tcgaCompare plot: x-axis labels getting clipped
I have a {maftools} plot:
whatever I try to save that completely I see the down part of plot is trimmed like

I have tried both save as pdf and save as image options which did not differ
Any suggestion please?
Plot comes by this code
tcgaCompare(maf = laml,cohortName = "Cohort",rm_zero = FALSE,decreasing = TRUE,logscale = TRUE,primarySite = TRUE,col = c("grey30", "red"),medianCol = "black")

One way is to increase the outer margins of plot.
par(oma=c(7,0,0,0))
tcgaCompare(maf = laml, primarySite = TRUE)


text:
Replacing some strings in R or terminal
I have a .txt file; I have 139 unique sample IDs in sampleID column
> head(mutations)
                                                                              sampleID chr     pos ref mut
4 LP2000107-DNA_A01_vs_LP2000102-DNA_A01__pv.1.6__rg.grch37_g1k__al.bwa_mem__.snp.pass   1  920552   G   A
5 LP2000107-DNA_A01_vs_LP2000102-DNA_A01__pv.1.6__rg.grch37_g1k__al.bwa_mem__.snp.pass   1 1350208   G   A
6 LP2000107-DNA_A01_vs_LP2000102-DNA_A01__pv.1.6__rg.grch37_g1k__al.bwa_mem__.snp.pass   1 1447367   T   G
7 LP2000107-DNA_A01_vs_LP2000102-DNA_A01__pv.1.6__rg.grch37_g1k__al.bwa_mem__.snp.pass   1 1909310   G   A
8 LP2000107-DNA_A01_vs_LP2000102-DNA_A01__pv.1.6__rg.grch37_g1k__al.bwa_mem__.snp.pass   1 2172675   G   C
9 LP2000107-DNA_A01_vs_LP2000102-DNA_A01__pv.1.6__rg.grch37_g1k__al.bwa_mem__.snp.pass   1 2846623   C   T
>

> dim(mutations)
[1] 3574407       5
>

I want to change name of sample from 1 to 139 in R or Linux terminal but I don't know how to do that, any help please?
Thanks

I believe that
uniqsampleID = 1:length(unique(mutations$sampleID))
names(uniqsampleID) = unique(mutations$sampleID)
mutations$sampleID = uniqsampleID[mutations$sampleID]

should work, assuming mutations is a data frame.

Explanation:
1:N creates a vector of length N, containing 1, 2, ..., N. So using 1:length(X) will do this for each item in X, so 1:length(unique(mutations$sampleID)) should give you a vector of 1:139.
Then names(uniqsampleID) = unique(mutations$sampleID) assigns to each of these numbers one of these sample IDs, which allows you to retrieve the sample number with uniqsampleID[mutations$sampleID].

Sorry that gives 3574407 sample IDs while I only have 139 samples for each I have 25 chromosome and a lot of mutations; I mean sampleID must contains maximum 139 unique numbers

Ah I see! I'll change it.

Yessssss worked thanks; Sorry in chr column I have the number of chromosome 1:23 plus to X and Y; How I can add chr prefix to this column? I mean instead of 1 I should have chr1

Technically should be another question, but paste0() can handle that. It would probably be mutations$chr = paste0("chr", mutations$chr)


text:
Getting weird plot with ggpubr package
I have a data frame for two groups and 9 variables
> head(aa)
  treatment variable     value
1       pre     Sig1 0.1173863
2       pre     Sig1 0.1888243
3       pre     Sig1 0.2191765
4       pre     Sig1 0.1277735
5       pre     Sig1 0.1711712
6       pre     Sig1 0.0000000
> class(aa)
[1] "data.frame"
> str(aa)
'data.frame':   1026 obs. of  3 variables:
 $ treatment: Factor w/ 2 levels "post","pre": 2 2 2 2 2 2 2 2 2 2 ...
 $ variable : Factor w/ 9 levels "Sig1","Sig2",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ value    : num  0.117 0.189 0.219 0.128 0.171 ...
> tail(aa)
     treatment variable value
1021      post     Sig9     0
1022      post     Sig9     0
1023      post     Sig9     0
1024      post     Sig9     0
1025      post     Sig9     0
1026      post     Sig9     0
>

I am using this code
library(ggpubr)
> p <- ggboxplot(aa, x = "variable", y = "value",color = "treatment")
> my_comparisons <- c("pre","post")
> p + stat_compare_means(comparisons = my_comparisons)+ stat_compare_means(label.y = 50)
Warning message:
Computation failed in `stat_signif()`:
missing value where TRUE/FALSE needed
>

But I am getting this weird plot

While something like this is my goal

I really don't know where I am doing wrong
I used ggplot where I got something like this but I need statistics that is why I have stuck on this

I then used this code
> DF <- aa %>% pivot_longer(.,cols = c(Sig1,Sig2,Sig3, Sig4,Sig5,Sig6,Sig7,Sig8,Sig9), names_to = "var", values_to = "val")
> View(DF)
> library(ggplot2)
> library(ggpubr)
> ggplot(DF, aes(x = treatment, y = val, fill = treatment))+
+     geom_boxplot()+
+     stat_compare_means(comparisons = c("pre","post")+
+     stat_compare_means()+
+     facet_grid(var~., scales = "free")
+

But I got another weird plot

But changing the scale I got this but gives one statistics

While I need one stat for each variable between pre and post groups
This solution gives to statistics
> p <- ggboxplot(aa, x = "variable", y = "value",color = "treatment",facet.by = "variable")
> p + stat_compare_means(comparisons = my_comparisons)+ stat_compare_means(label.y =2)

When I edited my code
> p <- ggboxplot(aa, x = "variable", y = "value",color = "treatment")
> p + stat_compare_means(aes(group = variable))
>

But again one stat

The problem is the scale used: For the plot you called "weird" (first from the top), the scale is 50 and for the "ggplot only" (third from the top) the scale is 1.
You should play with the stat_compare_means(label.y = 50) bit, you can try setting the label.y parameter to 1.5 or 2.
EDIT BECAUSE OF THE ADDITIONAL REMARKS ON THE QUESTION:
Once again (and for sure the last time), I had to manually create example data as OP ignores repeated reminders on how to make a reproducible example with R.
library(data.table)
library(ggpubr)

text <- "
treatment variable     value
1       pre     Sig1 0.1173863
2       pre     Sig1 0.1888243
3       pre     Sig1 0.2191765
4       pre     Sig1 0.1277735
5       pre     Sig1 0.1711712
6       pre     Sig1 0.0000000
1       post     Sig1 0.1173863
2       post     Sig1 0.1888243
3       post     Sig1 0.2191765
4       post     Sig1 0.1277735
5       post     Sig1 0.1711712
6       post     Sig1 0.0000000
1021      pre     Sig9     0
1022      pre     Sig9     0
1023      pre     Sig9     0
1024      pre     Sig9     0
1025      pre     Sig9     0
1026      pre     Sig9     0
1021      post     Sig9     1
1022      post     Sig9     1
1023      post     Sig9     1
1024      post     Sig9     1
1025      post     Sig9     1
1026      post     Sig9     1
"
my_data <-fread(text)

> my_data
      V1 treatment variable     value
 1:    1       pre     Sig1 0.1173863
 2:    2       pre     Sig1 0.1888243
 3:    3       pre     Sig1 0.2191765
 4:    4       pre     Sig1 0.1277735
 5:    5       pre     Sig1 0.1711712
 6:    6       pre     Sig1 0.0000000
 7:    1      post     Sig1 0.1173863
 8:    2      post     Sig1 0.1888243
 9:    3      post     Sig1 0.2191765
10:    4      post     Sig1 0.1277735
11:    5      post     Sig1 0.1711712
12:    6      post     Sig1 0.0000000
13: 1021       pre     Sig9 0.0000000
14: 1022       pre     Sig9 0.0000000
15: 1023       pre     Sig9 0.0000000
16: 1024       pre     Sig9 0.0000000
17: 1025       pre     Sig9 0.0000000
18: 1026       pre     Sig9 0.0000000
19: 1021      post     Sig9 1.0000000
20: 1022      post     Sig9 1.0000000
21: 1023      post     Sig9 1.0000000
22: 1024      post     Sig9 1.0000000
23: 1025      post     Sig9 1.0000000
24: 1026      post     Sig9 1.0000000
      V1 treatment variable     value

ggboxplot(my_data,
          x = "treatment",
          y = "value",
          color = "treatment",
          facet.by = "variable") +
  stat_compare_means(label = "p.format")

And this code generates:

Thank you. now seems better but gives one statistics for all, what should I do to get one statistics in comparing each variable between two group?

You will need to play with the comparisons parameter within the stat_compare_means() function.

Sorry how? I have two groups pre and post which I placed them in comparison

A workaround would be doing facet.by = "variable" within your ggboxplot() call.

Thank you but gives no stat at all

Why do you call stat_compare_means() twice? And what happens when you use stat_compare_means(aes(group = variable))?

I removed extra stat_compare_means(), and stat_compare_means(aes(group = variable)) gives one stat

Thank you so much, are the stats are Kruskal-wallis?


text:
Grouping a set of genes in a ggplot object
I have a boolean matrix of mutated or non mutated genes across some samples
> head(muts)
                  TP53 CDKN2A TSHZ3 ARID1A EPHA3 NOTCH1 ACVR2A ARID1B ARID2
LP6008336-DNA_G01    1      1     1      0     0      0      0      0     0
LP6008460-DNA_D01    1      1     0      1     0      0      0      0     0
LP6005500-DNA_D03    1      1     0      0     0      0      0      0     0
LP6008334-DNA_B02    1      0     1      0     0      0      0      0     0
LP6008334-DNA_A04    1      0     1      0     0      0      0      0     0
LP6005334-DNA_H01    1      0     0      0     1      0      0      1     0

By this code
pw_rt = data.frame(
  gene = c('PIK3CA', 'EGFR', 'ERBB2', 'PTEN'),
  pw = c('RTK'),
  stringsAsFactors = F
)

pw_cc = data.frame(
  gene = c('CDKN2A', 'CCNE1'),
  pw = c('Cell Cycle'),
  stringsAsFactors = F
)

pw = rbind(pw_rt, pw_cc)

pw_col = RColorBrewer::brewer.pal(n = 3, name = 'Set1')[1:2]
names(pw_col) = unique(pw$pw)

require(dplyr)

lmuts = reshape2::melt(muts) %>% as_tibble()
colnames(lmuts) = c('sample', 'gene', 'value')

pw_n = pw$pw
names(pw_n) = pw$gene

Np = nrow(muts)

lmuts %>%
  mutate(gene = paste(gene), PW = pw_n[gene]) %>%
  group_by(gene, PW) %>%
  summarise(N = sum(value == 1)) %>%
  ungroup() %>%
  mutate(
    PW = ifelse(is.na(PW), "None", PW),
    N = paste0(N, ' (', round(N/Np * 100, 1), '%)')
    ) %>%
  ggplot(aes(x = PW, y = gene, fill = PW)) +
  geom_tile() +
  geom_text(aes(label = N)) +
  theme_light() +
  theme(legend.position = 'bottom') +
  guides(fill = guide_legend('Pathway')) +
  scale_fill_manual(values = c(pw_col, `None` = 'gainsboro')) +
  scale_x_discrete(limits = c(names(pw_col), 'None')) +
  labs(
    x = "",
    y = 'Gene',
    title = "Occurrence of mutations"
  )

I am able to get proportion of altered pathway by these genes like

However I want to highlight oncogenes and tumour suppressors within this plot
For instance TP53 is a tumour suppressor and CCNE1 is an oncogene
I don't know how to manipulate the above code to highlight or partition genes to tumour suppressor or oncogene
Can you help?
EDITED
I have add extra column to lmuts like
> head(lmuts)
             sample gene value gene_class
1 LP6005690-DNA_H02 TP53     1        TSG
2 LP2000333-DNA_A01 TP53     1        TSG
3 LP6005409-DNA_D03 TP53     1        TSG
4 LP6008141-DNA_H02 TP53     1        TSG
5 LP6008336-DNA_E02 TP53     1        TSG
6 LP6008269-DNA_B06 TP53     1        TSG
>

but
> lmuts %>%
+     mutate(gene = paste(gene), PW = pw_n[gene]) %>%
+     group_by(gene, PW) %>%
+     summarise(N = sum(value == 1)) %>%
+     ungroup() %>%
+     mutate(
+         PW = ifelse(is.na(PW), "None", PW),
+         N = paste0(N, ' (', round(N/Np * 100, 1), '%)')
+     ) %>%
+     ggplot(aes(x = PW, y = gene, fill = PW)) +
+     geom_tile() +
+     geom_text(aes(label = N)) +
+     theme_light() +
+     theme(legend.position = 'bottom') +
+     guides(fill = guide_legend('Pathway')) +
+     scale_fill_manual(values = c(pw_col, `None` = 'gainsboro')) +
+     scale_x_discrete(limits = c(names(pw_col), 'None')) +
+     labs(
+         x = "",
+         y = 'Gene',
+         title = "Occurrence of mutations"
+     )+facet_wrap(~gene_class, ncol=1)
Error: At least one layer must contain all faceting variables: `gene_class`.
* Plot is missing `gene_class`
* Layer 1 is missing `gene_class`
* Layer 2 is missing `gene_class`
Run `rlang::last_error()` to see where the error occurred.

If you had the classification into oncogene/TSG for each gene in a column (let's call it gene_class) in lmuts, you could split both groups into facets within the plot by adding:
+ facet_wrap(~gene_class, ncol=1)

I can't reproduce your graph, because you don't make muts available, but it appears your problem is having x = PW in your ggplot call. Don't you want to that to be whatever your variable name for sample is (ie, the rownames in your but dataframe you have at the beginning).

Sorry I have edited my post, know you can reproduce my plot because muts is available now


text:
Batch correction in differential expression analysis
I have sent two sets (two batches in matter of sending for sequencing) of different samples (plasma) to small RNA-seq to Qiagen company
This is how my meta data look
> head(meta)
           batch group    sex
30678.009 batch2     C   Male
30678.010 batch2     C Female
30941.001 batch1     C Female
30941.002 batch1     C   Male
30941.003 batch1     C   Male
30941.004 batch1     C   Male
> tail(meta)
           batch group    sex
30941.016 batch1    LT Female
30941.017 batch1    LT   Male
30941.019 batch1    LT   Male
30941.022 batch1    LT Female
30941.023 batch1    LT   Male
30941.024 batch1    LT Female
>

I want to extract differentially expressed RNAs between my two groups without batch effect
I am doing this
Is this correct please?
dds <- DESeqDataSetFromMatrix(countData = count,
                              colData = meta,
                              design= ~ batch + group)

dds <- DESeq(dds, test="LRT", reduced=~batch)

res <- results(dds)
resultsNames(dds)

[1] "Intercept"              "batch_batch2_vs_batch1" "group_LT_vs_C"
>

If both batches contain replicates of both groups (LT and C) then the presented code is what you should do, exactly pointed out in the vignette for the LRT: http://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#likelihood-ratio-test

Thank you. Please, what do you mean by "If both batches contain replicates of both groups (LT and C)" ? I know each sample is one independent person and we have sent samples to Qiagen in two times

It means exactly like it reads. Both batches need to contain replicates (preferably in comparable numbers) of both groups. Use table(meta$batch, meta$group) to see the distribution contingency table.


text:
A two color bar plot in R
I have raw read counts of miRNAs in two conditions like below
> head(long1)
             miRNA Read count Condition
1:  hsa-miR-539-3p         33    living
2: hsa-miR-199b-5p        102    living
3: hsa-miR-146a-3p         27    living
4:   hsa-let-7a-3p        186    living
5:    hsa-miR-1275         54    living
6:  hsa-miR-29a-3p       6464    living
>

For each individual miRNAs, I want to show the difference of raw read counts across conditions by
p=ggplot(data = long1, aes(x=long1$miRNA, y=long1$`Read count`)) + geom_boxplot(aes(fill=Condition))
p <- p + guides(fill=guide_legend(title="Condition"))
p <- p + xlab("miRNAs") + ylab("Raw mapped reads") + ggtitle("s vs living with cancer")
p + facet_wrap( ~ miRNA, scales="free")+theme_bw() add_pval(p, pairs = list(c("s", "living")), test='wilcox.test')

Which returns this
[![enter image description here][1]][1]
This plot is not good
How I can a plot like below in which miRNAs are on Y axis and two colours shows the conditions like this

Using Google well is one of the best skillsets you can possibly learn in programming-allied fields. The style of plot you're looking for is a "grouped bar plot". I googled "r grouped bar plot" and found a nice link: https://www.r-graph-gallery.com/48-grouped-barplot-with-ggplot2.html
The key is to use position=dodge alongside color=Condition, x="Read count", y=miRNA )not sure ow you'd quote the name "Read count", might be better off renaming that column to readcount).


text:
Selecting a part of heatmap
I have a heatmap
I see a block of genes for a group of samples
How I can select a part of heatmap?
I mean is there any way to export a list of genes from rows of a heatmap in R?
I have tried which returns nothing
> library(ComplexHeatmap)

> ht = Heatmap(heat1)
> ht[1:35,1:5]

Why would you want to do that? Why not simply subset the data that generates the heatmap to the desired rows? What have you tried to do or searched to do this?

From a matrix of scaled values, how I can say which miRNAs I want before being clustered and visualized so that in the heatmap I can see a blocked of correlated miRNAs?

You can check/create clusters outside the heatmap via hclust. Keep in mind that the heatmap will change if you provided less genes/miRNAs or different samples. But probably can extract the data from ht object.

You should really read the manual before asking questions here.

You're using ComplexHeatmap, which is pretty easy to subset. Your subset operation creates a subset object of the heatmap, not of the data in the heatmap.
In case you don't want to run hclust outside of the CH, you can use row_order and column_order to re-order your data in the CH display order.
library(ComplexHeatmap)

ht = Heatmap(heat1)
r_o <- row_order(ht)
c_o <- column_order(ht)
heat1_reordered <- heat1[r_o, c_o]
heat1_reordered[1:35,1:5]

CH does say that clustering is oly "fixed" on drawing, so you may want to run this before picking the row and column orders:
pdf(NULL)
ht <- draw(ht)
dev.off()


text:
Looping over several files in bash
I have a bash script:
I am wondering how I can change this script to loop over a bunch of .vcf files and give output .txt with the name of corresponding .vcf
I tried changes done in similar script in  Changing this code in a way to work for my files but I just  messed everything up

If we save the script pasted in the main post as a sh file and we have some .vcf files in a folder, by this line we can iterate over vcf files to extract what mentioned in the script returning .txt output by the name of corresponding vcf
for file in *.vcf ; do bash indel_vcf_parasing.sh $file ; done

Can you provide some context to future readers about this answer

This is generally the method I use for scaling. Write a program to do 'one', then write wrappers, or use one-liners (as this answer shows) to do many.
You will need to be sure that successive runs do not step on each other (if using a single directory for all runs), or make a new directory for each run if necessary.

Please edit your answer to include the extra details so the answer is complete and a user doesn't have to read all the comments to understand.


text:
Calculating mean accross rows with repeated entries in R
I have this expression matrix
> head(dat[1:8,1:4])
  Gene.Symbol        p1        p2         p3
1       STAT1     7.6159     7.6784     7.4430
2       STAT1     7.4435     7.2850     7.1250
3       STAT1     7.8925     7.9269     7.4956
4       STAT1     7.5867     7.9285     7.6741
5       GAPDH     8.4861     8.3189     9.0913
6       GAPDH     9.0206     8.0813     9.2360
> dim(dat)
[1] 58324   276
>

As you are seeing, I have repeated genes, so how I take average over expression values for redundant genes finishing with one unique expression value for each gene for each sample in columns?

Using group_by from dplyr
You can use group_by function from dplyr to calculate the mean of each gene symbol in each samples:
library(dplyr)
dat %>%
  group_by(Gene.Symbol) %>%
  summarise_at(vars(contains("p")), funs(mean(., na.rm = TRUE)))

# A tibble: 2 x 4
  Gene.Symbol    p1    p2    p3
  <chr>       <dbl> <dbl> <dbl>
1 GAPDH        8.75  8.20  9.16
2 STAT1        7.63  7.70  7.43

Using data.table
Using data.table to speed up the calculation, you can do:
library(data.table)
setDT(dat)
dat[, lapply(.SD, mean), by = Gene.Symbol]

   Gene.Symbol      p1     p2       p3
1:       STAT1 7.63465 7.7047 7.434425
2:       GAPDH 8.75335 8.2001 9.163650
```


text:
Visualizing the difference among three groups
I have 3 data frames for three groups of patients and in each of them I have the number of mutation types like insertion, deletion, SNP and total of mutations for each patients. In each group I have different number of patients like below
> head(dat1)
                patient  DEL  INS   SNP   total
1:    LP6008337-DNA_H06  927  773 40756   42456
2:    LP6008334-DNA_D02 1049  799 31009   32857
> dim(dat1)
[1] 21  6
>

 head(dat2)
                Patient   DEL  INS   SNP total
1:    LP6008031-DNA_E01 13552 3374 62105   79031
2:    LP6005500-DNA_G01   539  500 43451   44490
> dim(dat2)
[1] 33  6
>

> head(dat3)
                Patient   DEL   INS    SNP    total
1:    LP6005935-DNA_F03 39168 16739  58095   114002
2:    LP6008269-DNA_D08   849   910 103501  105271

> dim(dat3)
[1] 106   6

I want to show if the number of each mutational category and total is different between these three groups by chi-square test
Actual question is if total number of mutations, SNP, DEL and INS are statistically significant among groups. I used pairwise t test but I afraid this test is not a suitable test for the distribution of data moreover I don't know how to visualize the p-value
Imagining the comparison of total number of mutations between three groups, this picture is a good example

Or in this plot they compare several features but only two groups

Can you help ?
Thanks

So make a chart with three groups. Where is the actual question?

Actual question is if total number of mutations, SNP, DEL, INS statistically significant among groups. I used pairwise t test but I afraid this is not a suitable test moreover I don't know how to visualize the p-value

Here a way to do it is to start first by creating a dataset containing your three different group, based on the few lines you display:
dat1$Condition = "dat1"
dat2$Condition = "dat2"
dat3$Condition = "dat3"
colnames(dat1)[1] = "Patient"
DF <- rbind(dat1,dat2,dat3)

             Patient   DEL   INS    SNP  total Condition
1: LP6008337-DNA_H06   927   773  40756  42456      dat1
2: LP6008334-DNA_D02  1049   799  31009  32857      dat1
3: LP6008031-DNA_E01 13552  3374  62105  79031      dat2
4: LP6005500-DNA_G01   539   500  43451  44490      dat2
5: LP6005935-DNA_F03 39168 16739  58095 114002      dat3
6: LP6008269-DNA_D08   849   910 103501 105271      dat3

Then, you can reshape DF in order to put all values DEL, INS, SNP and Total into a single column and create a categorical column. Doing this will be helpful for later producing the facets of your plot. To do that, you can use pivot_longer from tidyr package:
library(tidyr)
DF <- DF %>% pivot_longer(.,cols = c(DEL,INS,SNP, total), names_to = "var", values_to = "val")

# A tibble: 24 x 4
   Patient           Condition var     val
   <chr>             <chr>     <chr> <int>
 1 LP6008337-DNA_H06 dat1      DEL     927
 2 LP6008337-DNA_H06 dat1      INS     773
 3 LP6008337-DNA_H06 dat1      SNP   40756
 4 LP6008337-DNA_H06 dat1      total 42456
 5 LP6008334-DNA_D02 dat1      DEL    1049
 6 LP6008334-DNA_D02 dat1      INS     799
 7 LP6008334-DNA_D02 dat1      SNP   31009
 8 LP6008334-DNA_D02 dat1      total 32857
 9 LP6008031-DNA_E01 dat2      DEL   13552
10 LP6008031-DNA_E01 dat2      INS    3374
# … with 14 more rows

Now, you can plot your data and use the function stat_compare_means from ggpubr package to get the significance display on each tab. Based on the size of your group, you may want to adjust for the type of statistical test to be used (take a look at ?stat_compare_means to find the statistical test adequate to your dataset). You can pass various conditions to test as a list by using comparisons argument.
For facetting your graph, you can use facet_grid (orfacet_wrap, try the one you preferred). Altogether, you can get the plot by writing:
library(ggplot2)
library(ggpubr)
ggplot(DF, aes(x = Condition, y = val, fill = Condition))+
  geom_boxplot()+
  stat_compare_means(comparisons = list(c("dat1","dat2"), c("dat1","dat3"), c("dat2","dat3")))+
  stat_compare_means()+
  facet_grid(var~., scales = "free")

Hope it helps you to figure it out how to represent your data.

Thank you, you are amazingly helpful, how I can I do exactly the same procedure but for pairwise t-test instead of cruskal wallis? Because I found the pvalue completely different

Glad to be able to help. It's normal that Kruskal-Wallis returns different p values than pairwise.t.test because one is non-parametric and the other is parametric. I don't think you have the option to pass a parwise.t.test in stat_compare_means but you can try to add the argument method = "t.test" (e.g. stat_compare_means(method = "t.test")). But you should check that pairwise.t.test and t.test gives you the same output first. If not, you will have to add them manually using geom_signif or stat_pvalue_manual.

Thank you, this gives error . > ggplot(DF, aes(x = Response, y = val, fill = Response))+

geom_boxplot()+
stat_compare_means(comparisons = list(c("Responders","Non_res_naive"), c("Responders","Non_res_post"), c("Non_res_naive","Non_res_post")),method = "pairwise.t.test")+
stat_compare_means()+
facet_grid(var~., scales = "free")  + scale_y_continuous(trans = "log2")
Error in .method_info(method) :
Non-supported method specified. Allowed methods are one of: t.test, t.test, t.test, wilcox.test, wilcox.test, wilcox.test, anova, anova, kruskal.test, kruskal.test

That's because method = "pariwise.t.test does not exist in stat_compare_means, that's why, I advised you in my comment above to use method = "t.test"


text:
How to visualise the proportion of clonal genes to sub clonal genes
I have a data frame
Response_TO_CHEMOTHERAPY    Gene    CCF
Responders                  POLQ    1
Responders                 CDKN2A   1
Responders                  TP53    1
Responders              CCDC102B    0.906729075
Non-responders           CCDC102B   1
Non-responders           CCDC102B   1
Non-responders           CCDC102B   1
Non-responders             TSHZ3    1
Non-responders             PIK3CA   1
Responders                 SLIT2    1
Non-responders              SLIT2   1

I have two groups of patients; Responders to chemotherapy and non-responders. For each group I have calculated cancer cell faction (CCF) for a fixed set of genes so that an individual gene has several reads in each group
If CCF = 1 the gene clonal and if CCF < 1 gene is sub clonal, for a given gene for each group, I therefore would have the number of clonal and sub clonal for this gene, for example 2 clonal 7 sub clonal which would be 2/7 something like the below figure

This is another potential example which reveals what is in my mind

Which show the clonal frequency for each group
I don't know how to use dput() but this is my complete data
https://www.dropbox.com/s/ol7m4nlkuyuy85a/c1.txt?dl=0

Sorry why you are unvoting my question? Is not it clear enough? When I am asking a question that means I need to solve that not that I am posting a question for fun. I invest time to prepare input and to know what exactly I am looking for

I've used some dummy data for simplicity. It's fairly easy to do this with ggplot, but you have to format your matrix-like data for ggplot input.
eg :
# make data
set.seed(1)
ngenes <- 20
nsamps <- 10
my_dat <- matrix(rnorm(ngenes*nsamps), ncol = nsamps, nrow = ngenes)
colnames(my_dat) <- paste("samp", seq_len(nsamps))
rownames(my_dat) <- paste("gene", seq_len(ngenes))

library(ggplot2)
library(reshape2) # for melt

# format data for plotting
dat_toplot <- melt(my_dat)

#plot
ggplot(dat_toplot, aes(x = Var1, y = Var2, fill = value)) + geom_tile()

plots

If you need specifics on plot tweaking, google's your friend. There are plenty of tutorials on how to use ggplot for heatmaps. (eg. here, or here)
EDIT
I loaded your data into the dat object.
library(tidyverse)
dat <- as_tibble(dat)

dat_counts <- dat %>% group_by(Sample, Gene) %>% summarise(count_subclonal = sum(CCF<1), count_clonal = sum(CCF>=1))
dat_counts$ratio <- dat_counts$count_subclonal/dat_counts$count_clonal

# set infinite values (/0) to NA
dat_counts$ratio[is.infinite(dat_counts$ratio)] <- NA

ggplot(dat_counts, aes(x = Sample, y = Gene, fill = ratio)) + geom_tile() +
  geom_text(aes(label = paste(count_subclonal, "/", count_clonal)), col = 'white') # add ratio text

Outputs :

Thank you but this is not what I need, I need the proportion of clonal to sub clonal genes in each box not a heat map of cancer cell fraction

@Exhausted Ah, then your question is not on generating a heatmap. It is "how to compute the proportion of clonal to subclonal genes"

@Exhausted See my updated post


text:
Running a script in python
I am trying to run a script in python but I getting this error
loan-mac-13:vcf2maf fi1d18$ python
Python 2.7.10 (default, Feb 22 2019, 21:55:15)
[GCC 4.2.1 Compatible Apple LLVM 10.0.1 (clang-1001.0.37.14)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import os
>>> os.chdir("/Users/fi1d18/Downloads/vcf2maf-master/vcf2maf")
>>>
>>>
>>> os.getcwd()
'/Users/fi1d18/Downloads/vcf2maf-master/vcf2maf'
>>> vcf2maf.py -h
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
NameError: name 'vcf2maf' is not defined
>>> python vcf2maf.py -h
  File "<stdin>", line 1
    python vcf2maf.py -h
                 ^
SyntaxError: invalid syntax
>>>

Where I am doing wrong?
Thanks for any help
By your comments I run
Sorry says that Last login: Wed Apr  1 11:50:14 on ttys000
pc-133-235:~ fi1d18$ pwd
/Users/fi1d18
pc-133-235:~ fi1d18$
pc-133-235:~ fi1d18$ cd Downloads
pc-133-235:Downloads fi1d18$ python vcf2maf.py -h
Traceback (most recent call last):
  File "vcf2maf.py", line 19, in <module>
    import vcf
ImportError: No module named vcf
pc-133-235:Downloads fi1d18$

As I said, that's a completely separate issue. Please ask a new question about it. First, however, just copy/paste the error into your favorite search engine. That will give you the answer.

I am sorry to say, but this is a bioinformatics equivalent of "What is a cell" in biology. It appears to me that you really should go through some python tutorials to get over the basics before you proceed with your research. I heavily recommend https://software-carpentry.org/ as a starting point.

You %100 right @KamilSJaron, I am nothing in the first line of python

You are trying to run a script file in the python console. That doesn't make sense. The console is for running commands manually, interactively. A script is a collection of commands. You can't just run that by name in the console. If you want to use it, copy the contents of your script and try to run them in the console.
The way you run a python script file is simply:
python vcf2maf.py -h

Or, assuming you have made the file executable and it has the right shebang line, just:
./vcf2maf.py -h

The specific error you are getting is because foo.bar in Python means "run the bar method from the object foo. So vcf2maf.py means run py from vcf2maf. Since there is no object defined with the name vcf2maf, you get the name 'vcf2maf' is not defined.

Sorry says that Last login: Wed Apr  1 11:50:14 on ttys000
pc-133-235:~ fi1d18$ pwd
/Users/fi1d18
pc-133-235:~ fi1d18$
pc-133-235:~ fi1d18$ cd Downloads
pc-133-235:Downloads fi1d18$ python vcf2maf.py -h
Traceback (most recent call last):
File "vcf2maf.py", line 19, in
import vcf
ImportError: No module named vcf
pc-133-235:Downloads fi1d18$

@Angel I can't read that. Please [edit] your question to add extra information, don't put it in the comments. And tell me who/what "says" and when. I can't guess what commands you are running. If that's the output of running the script, then it is a completely different issue, you are missing the vcf module. So please ask a new question about that.