text:
Dodged bar chart
I want to create a dodged barchart with different cell cycles for each cell type within each condition.
dput(DF_count)
structure(list(orig.ident = c("treated", "treated", "treated",
"treated", "treated", "treated", "treated", "treated", "treated",
"treated", "treated", "treated", "treated", "treated", "treated",
"treated", "treated", "treated", "treated", "treated", "treated",
"treated", "treated", "treated", "treated", "treated", "treated",
"treated", "treated", "treated", "treated", "treated", "treated",
"untreated", "untreated", "untreated", "untreated", "untreated",
"untreated", "untreated", "untreated", "untreated", "untreated",
"untreated", "untreated", "untreated", "untreated", "untreated",
"untreated", "untreated", "untreated", "untreated", "untreated",
"untreated", "untreated", "untreated", "untreated", "untreated",
"untreated", "untreated", "untreated", "untreated", "untreated",
"untreated", "untreated", "untreated"), Phase = structure(c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L), .Label = c("G1", "G2M", "S"), class = "factor"), old.ident = structure(c(1L,
2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 1L, 2L, 3L, 4L, 5L,
6L, 7L, 8L, 9L, 10L, 11L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L,
10L, 11L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 1L, 2L,
3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 1L, 2L, 3L, 4L, 5L, 6L,
7L, 8L, 9L, 10L, 11L), .Label = c("Fibroblast", "T cell", "Macrophage",
"Stellate", "Acinar", "Endothelial", "Tumor", "B cell", "Mast cell",
"Ductal", "Islets of Langerhans"), class = "factor"), n = c(3257L,
666L, 560L, 446L, 91L, 464L, 192L, 506L, 162L, 207L, 90L, 603L,
608L, 238L, 98L, 22L, 110L, 215L, 95L, 63L, 97L, 10L, 2191L,
618L, 335L, 349L, 35L, 294L, 250L, 231L, 142L, 153L, 16L, 2910L,
1259L, 853L, 559L, 1134L, 348L, 306L, 133L, 116L, 12L, 192L,
706L, 1186L, 374L, 112L, 616L, 106L, 277L, 32L, 47L, 6L, 27L,
1715L, 1312L, 575L, 390L, 1036L, 250L, 252L, 91L, 66L, 8L, 68L
), Freq = c(53.8258139150554, 35.2008456659619, 49.4263018534863,
49.9440089585666, 61.4864864864865, 53.4562211981567, 29.2237442922374,
60.8173076923077, 44.141689373297, 45.2954048140044, 77.5862068965517,
9.96529499256321, 32.1353065539112, 21.0061782877317, 10.9742441209406,
14.8648648648649, 12.6728110599078, 32.7245053272451, 11.4182692307692,
17.1662125340599, 21.2253829321663, 8.62068965517241, 36.2088910923814,
32.6638477801269, 29.567519858782, 39.0817469204927, 23.6486486486486,
33.8709677419355, 38.0517503805175, 27.7644230769231, 38.6920980926431,
33.4792122538293, 13.7931034482759, 54.5863815419246, 33.5107798775619,
47.3362930077692, 52.6861451460886, 40.7035175879397, 49.4318181818182,
36.6467065868263, 51.953125, 50.6550218340611, 46.1538461538462,
66.8989547038328, 13.2432939410992, 31.5677402182592, 20.7547169811321,
10.5560791705938, 22.1105527638191, 15.0568181818182, 33.1736526946108,
12.5, 20.5240174672489, 23.0769230769231, 9.40766550522648, 32.1703245169762,
34.9214799041789, 31.9089900110988, 36.7577756833176, 37.1859296482412,
35.5113636363636, 30.1796407185629, 35.546875, 28.82096069869,
30.7692307692308, 23.6933797909408)), class = c("grouped_df",
"tbl_df", "tbl", "data.frame"), row.names = c(NA, -66L), groups = structure(list(
    orig.ident = c("treated", "treated", "treated", "treated",
    "treated", "treated", "treated", "treated", "treated", "treated",
    "treated", "untreated", "untreated", "untreated", "untreated",
    "untreated", "untreated", "untreated", "untreated", "untreated",
    "untreated", "untreated"), old.ident = structure(c(1L, 2L,
    3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 1L, 2L, 3L, 4L, 5L,
    6L, 7L, 8L, 9L, 10L, 11L), .Label = c("Fibroblast", "T cell",
    "Macrophage", "Stellate", "Acinar", "Endothelial", "Tumor",
    "B cell", "Mast cell", "Ductal", "Islets of Langerhans"), class = "factor"),
    .rows = list(c(1L, 12L, 23L), c(2L, 13L, 24L), c(3L, 14L,
    25L), c(4L, 15L, 26L), c(5L, 16L, 27L), c(6L, 17L, 28L),
        c(7L, 18L, 29L), c(8L, 19L, 30L), c(9L, 20L, 31L), c(10L,
        21L, 32L), c(11L, 22L, 33L), c(34L, 45L, 56L), c(35L,
        46L, 57L), c(36L, 47L, 58L), c(37L, 48L, 59L), c(38L,
        49L, 60L), c(39L, 50L, 61L), c(40L, 51L, 62L), c(41L,
        52L, 63L), c(42L, 53L, 64L), c(43L, 54L, 65L), c(44L,
        55L, 66L))), row.names = c(NA, -22L), class = c("tbl_df",
"tbl", "data.frame"), .drop = TRUE))

Used this code to generate a stacked bar graph:
ggplot(DF_Count, aes(x = old.ident, y = Freq, fill = Phase))+
  geom_col()+
  geom_text(aes(label = paste(round(Freq, 2),"%")),position = position_stack(vjust = 0.5))+
  facet_wrap(~orig.ident)

But want the 2 conditions to be next to each other on the same graph. How do I use this dataset and create a dodged bar chart to have both conditions together?

Here is a rather straight-forward work-around for your problem:
ggplot(DF_Count, aes(x = orig.ident, y = Freq, fill = Phase))+
  geom_col()+
  geom_text(aes(label = paste(round(Freq, 2),"%")),position = position_stack(vjust = 0.5))+
  facet_wrap(~old.ident)

It does not give a single plot or a dodged bar chart representing the two conditions you would like to visualize. For an exact solution to your problem, see this (actually the answer below the linked answer is analogous to what I have suggested above).


text:
Bar plot with error bars
I am trying to plot a box plot/stacked bar graph with error bars.
output of dput(DF)
`structure(list(Cell_Type = structure(c(1L, 2L, 3L, 4L, 5L, 6L,
7L, 1L, 2L, 3L, 4L, 5L, 6L, 7L), .Label = c("Fibroblast", "T cell",
"Macrophage", "Tumor", "Islets of Langerhans", "Endothelial",
"B cell"), class = "factor"), Condition = structure(c(1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("treated",
"untreated"), class = "factor"), Freq = c(6051L, 1892L, 1133L,
657L, 116L, 868L, 832L, 5331L, 3757L, 1802L, 835L, 287L, 704L,
256L), Percent = c(52.3941466793662, 16.3823707680319, 9.8103731924842,
5.68880422547407, 1.00441596675037, 7.51580223395965, 7.20408693393367,
41.0962072155412, 28.9623805118717, 13.8914585260561, 6.43694110391613,
2.21245760098674, 5.42707369719396, 1.97348134443417)), class = c("grouped_df",
"tbl_df", "tbl", "data.frame"), row.names = c(NA, -14L), groups = structure(list(
    Condition = structure(1:2, .Label = c("treated", "untreated"
    ), class = "factor"), .rows = list(1:7, 8:14)), row.names = c(NA,
-2L), class = c("tbl_df", "tbl", "data.frame"), .drop = TRUE))`

Gives me the following dataset.
       Var1      Var2 Freq
1   Fibroblast   treated 6051
2       T cell   treated 1892
3   Macrophage   treated 1133
4     Stellate   treated  893
5       Acinar   treated  148
6  Endothelial   treated  868
7   Fibroblast untreated 5331
8       T cell untreated 3757
9   Macrophage untreated 1802
10    Stellate untreated 1061
11      Acinar untreated 2786
12 Endothelial untreated  704

library(dplyr)
DF <- df %>% rename(Cell_Type = Var1, Condition = Var2) %>%
  group_by(Condition) %>%
  mutate(Percent = Freq / sum(Freq)*100)

# A tibble: 12 x 4
# Groups:   Condition [2]
   Cell_Type   Condition  Freq Percent
   <fct>       <fct>     <int>   <dbl>
 1 Fibroblast  treated    6051   55.1
 2 T cell      treated    1892   17.2
 3 Macrophage  treated    1133   10.3
 4 Stellate    treated     893    8.13
 5 Acinar      treated     148    1.35
 6 Endothelial treated     868    7.90
 7 Fibroblast  untreated  5331   34.5
 8 T cell      untreated  3757   24.3
 9 Macrophage  untreated  1802   11.7
10 Stellate    untreated  1061    6.87
11 Acinar      untreated  2786   18.0
12 Endothelial untreated   704    4.56

Calculated SD/Mean with this:
my_sum_1 <- DF %>%
    group_by(Cell_Type) %>%
    summarise(
        n=n(),
        mean=mean(Freq),
        sd=sd(Freq)
    ) %>%
    mutate( se=sd/sqrt(n))  %>%
    mutate( ic=se * qt((1-0.05)/2 + .5, n-1))

Then plotted with this:
ggplot(my_sum_1) +
    geom_bar( aes(x=Cell_Type, y=mean), stat="identity", fill="forestgreen", alpha=0.5) +
    geom_errorbar( aes(x=Cell_Type, ymin=mean-ic, ymax=mean+ic), width=0.4, colour="orange", alpha=0.9, size=1.5) +
    ggtitle("Cell Types")

I want to plot based on condition. Not sure how I would do that. I tried the fill command but doesn't seem to work as the object is not found, I can't seem to add it back to the data set. How do I include the condition in the bar chart?
When I'm calculating sd, I think it's combining the cell type value into one, rather than taking it individually. How do I change that?
dput(head(all@meta.data))
structure(list(orig.ident = c("PDAC", "PDAC", "PDAC", "PDAC",
"PDAC", "PDAC"), nCount_RNA = c(7945, 7616, 7849, 5499, 853,
1039), nFeature_RNA = c(2497L, 2272L, 2303L, 2229L, 509L, 588L
), percent.mt = c(2.63058527375708, 3.7421218487395, 4.9433048796025,
0.490998363338789, 0.234466588511137, 0.192492781520693), RNA_snn_res.0.5 = structure(c(5L,
5L, 5L, 5L, 1L, 6L), .Label = c("0", "1", "2", "3", "4", "5",
"6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16",
"17", "18", "19"), class = "factor"), seurat_clusters = structure(c(5L,
5L, 5L, 5L, 1L, 6L), .Label = c("0", "1", "2", "3", "4", "5",
"6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16",
"17", "18", "19"), class = "factor"), S.Score = c(0.0171539298241524,
-0.00599651468836135, 0.011058574403656, -0.0424413740104151,
0.03549521163095, -0.0243347616753686), G2M.Score = c(-0.0591387992008088,
0.0253275642795205, -0.0512402869839816, -0.0239076248512967,
-0.0263164126515597, -0.0339086517371782), Phase = structure(c(3L,
2L, 3L, 1L, 3L, 1L), .Label = c("G1", "G2M", "S"), class = "factor"),
    old.ident = structure(c(1L, 1L, 1L, 1L, 1L, 3L), .Label = c("Fibroblast",
    "T Cell", "Endothelial", "Tumor", "Stellate", "Macrophage",
    "B Cell", "Mast Cell", "Acinar", "Endocrine", "Exocrine"), class = "factor")), row.names = c("PDAC_AAACCCAGTCGGCTAC-1",
"PDAC_AAACGAAGTCCAGGTC-1", "PDAC_AAAGAACAGCAAGTGC-1", "PDAC_AAAGAACAGGATGAGA-1",
"PDAC_AAAGTGACATCCGTTC-1", "PDAC_AACAACCAGGAAGTAG-1"), class = "data.frame")

Your issue is coming by the fact that when you calculate mean, se, ... you removed the Condition column, here is your my_sum_1 results:
# A tibble: 7 x 6
  Cell_Type                n  mean    sd    se     ic
  <fct>                <int> <dbl> <dbl> <dbl>  <dbl>
1 Fibroblast               2 5691   509. 360    4574.
2 T cell                   2 2824. 1319. 932.  11849.
3 Macrophage               2 1468.  473. 334.   4250.
4 Tumor                    2  746   126.  89.0  1131.
5 Islets of Langerhans     2  202.  121.  85.5  1086.
6 Endothelial              2  786   116.  82    1042.
7 B cell                   2  544   407. 288    3659.

To calculate the mean, sd, ... from each cell types in each conditions, you need to group by Cell_Type and Condition by writing:
my_sum_1 <- DF %>%
  group_by(Cell_Type, Condition) %>%
  summarise(
    n=n(),
    mean=mean(Freq),
    sd=sd(Freq)
  ) %>%
  mutate( se=sd/sqrt(n))  %>%
  mutate( ic=se * qt((1-0.05)/2 + .5, n-1))

But now, your table looks like:
# A tibble: 14 x 7
# Groups:   Cell_Type [7]
   Cell_Type            Condition     n  mean    sd    se    ic
   <fct>                <fct>     <int> <dbl> <dbl> <dbl> <dbl>
 1 Fibroblast           treated       1  6051    NA    NA    NA
 2 Fibroblast           untreated     1  5331    NA    NA    NA
 3 T cell               treated       1  1892    NA    NA    NA
 4 T cell               untreated     1  3757    NA    NA    NA
 5 Macrophage           treated       1  1133    NA    NA    NA
 6 Macrophage           untreated     1  1802    NA    NA    NA
 7 Tumor                treated       1   657    NA    NA    NA
 8 Tumor                untreated     1   835    NA    NA    NA
 9 Islets of Langerhans treated       1   116    NA    NA    NA
10 Islets of Langerhans untreated     1   287    NA    NA    NA
11 Endothelial          treated       1   868    NA    NA    NA
12 Endothelial          untreated     1   704    NA    NA    NA
13 B cell               treated       1   832    NA    NA    NA
14 B cell               untreated     1   256    NA    NA    NA

Because in your initial dataset, you only have one value for each cell type in each condition, so you can calculate sd, sem. You need to first merge your replicates into your dataframe in order to be able to calculate your sem.
Then, you will be able to do:
library(ggplot2)
ggplot(my_sum_1, aes(x=Cell_Type, y=mean, fill = Condition)) +
  geom_col(alpha=0.5) +
  geom_errorbar( aes(ymin=mean-ic, ymax=mean+ic), width=0.4, colour="orange", alpha=0.9, size=1.5) +
  ggtitle("Cell Types")

Hope it helps you to figure it out what is your issue with your current data.

How do I merge a dataframe? These aren't replicates though. For each cell type, I have two conditions.

If you made replicates of your experiment, you should have multiple tables similar to  DF (one for each replicate) with the count of each cell types per conditions. If so, you can merge all of them (take a look for function ?merge() or ?full_join() from dplyr (https://dplyr.tidyverse.org/reference/join.html)). If all of your tables have the same format and colnames, you can consider using rbind (and adding first a condition column to each table corresponding to the replicate number).

I uploaded the question with the dataset before I calculated the frequencies and percentages.

I don't think it changes anything. This table does not represent replicates right ? Do you have replicates for this experiment ? What I mean by replicates is that you run a scRNAseq several times on your biological samples, and you identify the percentage of cell types in each of this runs. How do you design your experiments ? Can you describe it a little bit ? (no need to know exact parameters, just what you are trying to do and what materials did you use at the beginning)

I run one sample on the 10x chromium box. This data set contains samples from treated and untreated patients. 15 treated and 16 untreated samples. I run all the raw files through Seurat.

I uploaded the metadata from Seurat.

source: Stack Exchange

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

sort: votes


text:
Making a volcano plot look less cluttered
I have generated a volcano plot with a differential expression file.
Code for inputing file:
library(tidyr)
final_tumor <- df %>% pivot_longer(., -c(Feature.ID, Feature.Name), names_to = c("set",".value"), names_pattern = "(.+)_(.+)")

# A tibble: 80 x 6
Feature.ID Feature.Name set       Mean.Counts Log2.fold.change Adjusted.p.value
<fct>      <fct>        <chr>           <dbl>            <dbl>            <dbl>
  1 a          A            Cluster.1    0.000961            0.292           1
2 a          A            Cluster.2    0.000902            0.793           1
3 a          A            Cluster.3    0.00181             1.46            0.758
4 a          A            Cluster.4    0.000642            0.269           1
5 b          B            Cluster.1    0.000320            1.95            0.910
6 b          B            Cluster.2    0.00180             4.77            0.154
7 b          B            Cluster.3    0                   2.19            1
8 b          B            Cluster.4    0                   1.66            1
9 c          C            Cluster.1    0.00128            -2.01            0.0467
10 c          C            Cluster.2    0.00632             0.352           1
# … with 70 more rows

Here the output of head(final_tumor)
# A tibble: 6 x 6
Feature.ID        Feature.Name set       Mean.Counts Log2.fold.change Adjusted.p.value
<fct>             <fct>        <chr>           <dbl>            <dbl>            <dbl>
  1 ENSG00000227232.5 WASH7P       Cluster.1     0                   1.50            1
2 ENSG00000227232.5 WASH7P       Cluster.2     0                   1.73            1
3 ENSG00000227232.5 WASH7P       Cluster.3     0                   1.77            1
4 ENSG00000227232.5 WASH7P       Cluster.4     0.00114             4.30            0.293
5 ENSG00000227232.5 WASH7P       Cluster.5     0                   2.15            1
6 ENSG00000227232.5 WASH7P       Cluster.6     0                   1.22            1

And here is the output of tail(final_tumor):
# A tibble: 6 x 6
Feature.ID        Feature.Name set        Mean.Counts Log2.fold.change Adjusted.p.value
<fct>             <fct>        <chr>            <dbl>            <dbl>            <dbl>
  1 ENSG00000210196.2 MT-TP        Cluster.6       0.0699          -0.202           0.790
2 ENSG00000210196.2 MT-TP        Cluster.7       0.0801           0.0386          1
3 ENSG00000210196.2 MT-TP        Cluster.8       0.0711           0.0875          1
4 ENSG00000210196.2 MT-TP        Cluster.9       0.0152          -2.31            0.00127
5 ENSG00000210196.2 MT-TP        Cluster.10      0.0147          -2.30            0.00612
6 ENSG00000210196.2 MT-TP        Cluster.11      0.122            0.762           1

Here is the code for generating volcano plot:
library(ggplot2)
library(ggrepel)
ggplot(final_tumor, aes(x = Log2.fold.change,y = -log10(Adjusted.p.value), label = Feature.Name))+
  geom_point()+
  geom_text_repel(data = subset(final_tumor, Adjusted.p.value < 0.05),
                  aes(label = Feature.Name))+
  geom_text_repel(data = subset(final_tumor, Feature.Name == "Casp14"),
                  aes(label = Feature.Name), color = "red")

Generating the ggplot with all the data was taking forever, it was suggested to select random rows of the dataframe:
rep_final = final_tumor[sample(1:nrow(final_tumor), size = 1000), ]

Then I did:
ggplot(rep_final, aes(x = Log2.fold.change,y = -log10(Adjusted.p.value), label = Feature.Name))+
  geom_point()+
  geom_text_repel(data = subset(rep_final, Adjusted.p.value < 0.05), aes(label = Feature.Name))+
  geom_text_repel(data = subset(final_tumor, Feature.Name == "Casp14"), aes(label = Feature.Name), color = "red")

Everything is very tight and cluttered. I want to show the top differentially expressed genes. How do I go about simplifying the plot?
Is there a way to do this on EnhancedVolcano, pulling out a certain gene?
EnhancedVolcano(final_tumor, lab = as.character(final_tumor$FeatureName), x = 'Log2.fold.change', y = 'Adjusted.p.value', xlim = c(-8,8), title = 'Tumor', pCutoff = 10e-5, FCcutoff = 1.5, pointSize = 3.0, labSize = 3.0)

The element that is clustering your plot the most is the ggrepel labels, so that's what you need to work on.
To get to a readable number of labels, I would not remove random genes: you risk losing important info, be it genes that could be of interest to you or some of the most regulated. Incidentally that seems to happen on your example: you try plotting in red the label for Casp14, but it does not show on the plot, probably because the gene was removed in your sampling.
What you can do however is, as dc37 suggests, increase the thresholds for the labels (you can probably play around with the thresholds to see when it does become readable). This should also help with the plotting: you mention that generating the plot was taking a long time, but with the number of genes you are looking at (~10k-30k?), the geom_points() should not be a real problem (maybe a few seconds), however the ggrepel is the most computationnally expensive so that's what you should try to reduce.
This should get you started:
ggplot(final_tumor, aes(x = Log2.fold.change,y = -log10(Adjusted.p.value), label = Feature.Name))+
  geom_point()+
  geom_text_repel(data = subset(final_tumor, Adjusted.p.value < 0.001 & abs(Log2.fold.change) > 1 ), aes(label = Feature.Name))+
  geom_text_repel(data = subset(final_tumor, Feature.Name == "Casp14"), aes(label = Feature.Name), color = "red")

Notice everything plots the "final_tumor" object and not the "rep_final". Feel free to be even more stringent with the cutoffs if it's still taking too long to plot.
A bit outside of the score of your question, but I also wanted to mention that you should never be trying to label too many genes on a volcano plot: this type of graph is useful to visualize your data's p-value/FC distribution and highlight a subset of genes in this context (most regulated, genes regulated in another experiment, pathway of interest...), but if you label too many of them it will never be readable.
If you really want to be able to indentify all the genes on your graph, I would try an interactive visualization (like with plotly). I'm not going to detail the code here but you could look in that direction! (However, as with ggrepel, you might need to remove some genes so that it does not "take forever").
Have fun :)

It's still taking quite a bit of time to plot.

@mmpp how long?

It's still going.

Any other suggestions?

Aside from the ggrepel labels, another reason your volcano plot looks cluttered is the presence of the 4 very highly significant genes on the top right, which lead to the rest of the genes being squeezed into the x axis. For visualization purposes, you could try applying an inverse hyperbolic sine transformation (asinh() in R) to both your x/y axes, which should alleviate this squeezing while not distorting too much the overall volcano shape. Just be sure to also apply this transformation to any cutoff lines you are including in the plot.

As we discussed in comments in your previous post: Pulling out a certain gene in a volcano plot, in order to have less labeled genes on your volcano plot, you can specify more stringent cutoff.
For example, you can decide to display only genes with an absolute fold change superior to 2 by selecting the appropriate subset of your data:
ggplot(rep_final, aes(x = Log2.fold.change,y = -log10(Adjusted.p.value), label = Feature.Name))+
  geom_point()+
  geom_text_repel(data = subset(rep_final, abs(Log2.fold.change) > 2), aes(label = Feature.Name))+
  geom_text_repel(data = subset(final_tumor, Feature.Name == "Casp14"), aes(label = Feature.Name),
                  color = "red", size = 5)

Alternatively, you can maybe decide that a p value inferior to 0.05 is not enough stringent and decide to labeled only genes with an adjusted p value inferior to 0.00001 by doing:
ggplot(rep_final, aes(x = Log2.fold.change,y = -log10(Adjusted.p.value), label = Feature.Name))+
  geom_point()+
  geom_text_repel(data = subset(rep_final, Adjusted.p.value < 0.00001), aes(label = Feature.Name))+
  geom_text_repel(data = subset(final_tumor, Feature.Name == "Casp14"), aes(label = Feature.Name),
                  color = "red", size = 5)

Or, you can want to do a combination of both criteria by doing:
ggplot(rep_final, aes(x = Log2.fold.change,y = -log10(Adjusted.p.value), label = Feature.Name))+
  geom_point()+
  geom_text_repel(data = subset(rep_final, abs(Log2.fold.change) > 2 & Adjusted.p.value < 0.00001), aes(label = Feature.Name))+
  geom_text_repel(data = subset(final_tumor, Feature.Name == "Casp14"), aes(label = Feature.Name),
                  color = "red", size = 5)

At the end, I think it depends of what you are trying to demonstrate by presenting this volcano plot. Do you want to highlight some particular genes? Or do you want to display genes that are the most differentially expressed? If so, how many?
All of this depends of the scientific question underlying your experiments and why you came to make this volcano plot.
Hope this couple of example helps you to figure it out how to customize your volcano plot.

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


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 calculate coverage of each base of mapped reads?
I have a RNA aligned file without any header info.
I need to calculate the coverage of each mapped base for per chromosome. Is it possible to do it with perl ?
NS500377:28:H16A7BGXX:2:21205:16038:7658        99      chr2   3000482 255     46M     =       3000502 56      GTCCTACAGTGTGCATTTCTCGTTTTCCACGTTTCT    AAAA)FFFFFFFFFAFFFFFFFFFFF.F.FFFFFAF    XA:i:2  MD:Z:21T4T7T1   NM:i:3
NS500377:28:H16A7BGXX:3:11503:24953:10427       99      chr2   3000482 255     50M     =       3000513 71      GTCCTACAGTGTGCATTTCTCATTTTTCAAGTTTTTCAGT        AAAAAFFFFFFFFFFFFAFFFFFF7AAFFFFFFFFFFFFF
        XA:i:1  MD:Z:21T7C10    NM:i:2
NS500377:28:H16A7BGXX:3:12404:23929:9392        163     chr2   3000482 255     49M     =       3000483 40      GTCCTACAGTGTGCATTTCCCTTTTTTCACGTTTTTTAG AAFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF XA:i:1  MD:Z:19T16C2    NM:i:2
NS500377:28:H16A7BGXX:3:12404:23929:9392        83      chr2   3000483 255     49M     =       3000482 -40     TCCTACAGTGTGCATTTCCCTTTTTTCACGTTTTTTAGT AFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF XA:i:2  MD:Z:18T16C3    NM:i:2
NS500377:28:H16A7BGXX:3:21606:15675:1639        83      chr2   3000483 255     50M     =       3000480 -43     TCTTACAGTGTGCATTTCTCTTTTTTCACGTTTTTCAGTG        AFFFFFFFFFAFFFFFFFFFFFFFFFFFFFFFFFFAAAAA
        XA:i:0  MD:Z:2C37       NM:i:1
NS500377:28:H16A7BGXX:4:12605:18505:11088       83      chr2   3000486 255     50M     =       3000481 -45     TACAATGTGCATTTCTCTTTTTTCACGTTTTTCAGTATTT        FFFFFFFAFFF)FFFFFFFFF<FFFF<AFFFFFFFAA<AA
        XA:i:1  MD:Z:4G31G3     NM:i:2
NS500377:28:H16A7BGXX:2:13303:20634:12085       99      chr2   3000487 255     49M     =       3000521 73      ACAGTGTACATTTCTCATTTTTCACGTTTTTCAGTGTTA AAAAAFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF. XA:i:2  MD:Z:7G8T21T0   NM:i:3

You are looking for bedtools genomecov, specifically -d option. See https://bedtools.readthedocs.io/en/latest/content/tools/genomecov.html , for example:
bedtools genomecov -i A.bed -g my.genome -d

I recommend not reinventing the wheel: use bedtools instead of reimplementing the same functionality in Perl.

thanks, but bedtools doesn't work in this case as there is no header info.

Can you add it from another bam file using samtools? Or add a fake header?

i did add a header from a different sample to calculate the coverage... just wondering if we can do it using awk.


text:
Ancestry of the coronavirus 2019-nCov, WuHan city, China
In one of the answers to another question about the corona virus a link was given to this phylogenetic analysis of the virus.
Can somebody assist a non-bio type here? It seems to show that the current corona virus split from a virus in bats. And that the same ancestral virus is also ancestral to several other viruses, some in humans and some in bats.
Is my understanding basically correct?

I agree with your analysis, albeit this is a parsimonious analysis, but there is more going on here.

Overview
The central focus of the tree is to highlight the key biological concern of the new coronavirus, 2019-nCov. The key concern is the genetic similarities to SARS epidemic, and relates to the SARS receptor.
SARS background SARS is endemic in bats (your BioRxiv tree partly shows that and this tree definitely shows it) and in the 2002 epidemic infected civet cats which then infected humans. More importantly it is then transmitted from one human to another. This is a consequence of the SARS receptor been able to exploit the cellular receptors in the respiratory tract of bats, civet cats and humans and use this tissue as a site of replication. The concern is exacerbated because the tree here shows SARS independently infected humans on two separate occasions, suggesting the cross-over has an underlying genetic basis, which is part of the "SARS lineage".
BioRxiv 2019-nCov tree Your tree shows the 2019-nCov has a recent common ancestor with SARS in comparison to the rest of the betacoronaviruses. This therefore provides some circumstantial evidence the receptor mechanism and ability to frequently crossover from bats and ultimately transmit human to human could be shared with SARS. At an evolutionary level (which is what ANY tree is supposed to reconstruct) it raises the question which common ancestor of SARS did the "SARS receptor" originate. I'd need to draw a diagram to better demonstrate this point, but I hope you get the idea.
The "influenza receptor" The analogy is with influenza virus and the entire epidemiology of influenza hinges around is sialic acid receptor of influenza and its ability to bind on to the cellular receptors in the upper/lower respiratory tract of birds, pigs and humans. The hypothesised mechanism is called the mixing vessel theory and is the classic epidemiological understanding of how new pathogenic influenza pandemics occur. If you replace "birds and pigs" with e.g. "bats, civet cats and humans" you get the idea why this could be scary. We don't the intemediate host of 2019-nCov ... but I speculate there must be one, unless eating Cov-infected bats is common in China.
Technical details of the BioRxiv tree The tree is a nicely diverse selection of the beta-coronaviruses. The authors have rooted the tree using the outgroups of delacoronaviruses and gammacoronaviruses, so it is a good robust selection of outgroups which can be used to correctly identify the direct of evolution of betacoronavirus divergence. In tree theory (its formal name is phylogenetics theory), extensive rooting is good and minimises artefacts.
One of the earliest members of the betacoronaviruses to diverge is MERS (Middle Eastern respiratory syndrome), which form a single "clade" (all viruses share a unique common ancestor) and these represent around 50% amino acid divergence from SARS/2019-nCov. The selection of MERS in the tree from both camels and humans looks good. The selection of other betacoronaviruses looks great, I wasn't aware of the "ruminant clade"  at all involving buffalo, cow etc .. infections and there was an associated human infection. There are loads of bat isolates of the betacoronaviruses throughout the tree, but we eventually arrive at the "SARS clade". The authors show that 2019-nCov is an outgroup to the SARS clade and shows a close relationship to one (BioRxiv tree) or two (this tree) bat isolates. Looking at the precise SARS clade (better shown in this tree) we can see loads of bat virus associating with SARS lineages. We therefore assume the reservoir to SARS and likely 2019-nCov is bat, moreover that the single ancestor to both viruses was a bat (it is called a parsimonious hypothesis).
The one thing the BioRXiv tree omits, for example with regards this tree, is the diversity of SARS and in particularly the two independent origins of SARS, which is a weakness of their analysis, particular if this fed into downstream analysis. It is not to say the authors were wrong, but it was uncool.
Word of caution my understanding is the divergence between SARS and 2010-nCov is around 15% and this is a quite a large amount of genetic divergence, even if both viruses share a most recent common-ancestor, form part of the same clade and receptor peptide motifs. Nevertheless there is sufficient amino acid divergence to generate notable differences in epidemiology, clinical symptom and transmission.
MERS, SARS and 2019-nCov mortality rates It is worth noting that within the betacoronaviruses the mortality rates between difference clinically important viruses is very different. MERS has a mortality rate of 40-60%, SARS is around 10% but 2019-nCov is <2.5%. The 2019-nCov mortality rate is still important given it has infected more and could potentially infect alot more people than SARS.

Eating bats is common in the area of Wuhan.  WARNING! Potentially upsetting images. Google search for "bat soup."

Thanks @sockpuppet, :-( at least it they are boiled and the virus denatured :-(

@MichaelG - I presume your 2019-nCov mortality rate is deaths / cases. At this early stage of the outbreak, while the virus has not run it's course in the large majority of patients, I think a better estimate of the mortality rate is deaths / (deaths + recovered). Current data puts that at about 33%. Given the daily fluctuations in that quotient, I'd wager that in 6 months the mortality rate will be widely considered to be between 20 to 45%. Am I wrong about the calculation? Why the 2.5%?

@Mick I agree it is too early to guesstimate the mortality rate. My personal problem with the 2 to 2.5% mortality estimate is that it only applies to hospitalized cases and mild infections are not recorded. This sampling bias is very common in acute viral infections eg. flaviviruses, and would otherwise depress the mortality rate further. Your approach is ok for Ebolaviris, but I urge against high rate mortality theoretical approaches.

i heard that sars2 can not to live in bats because their hight temperature would damnege virus... so this virus has one story of years on something not a bat, not human or virus is ogm

Do you have a reference for this study @RosLuP? That would be interesting, at present we would assume it is a horseshoe bat and direct transmission to humans.

i remember someone show one protein max temperature and said that can not live in one super hot bat, or something as that... nothing scientific... i not remember where....

I came back to review this whole thing just now. My comment about "bat soup" is incorrect. The images you may see if you Google for "bat soup" are not from China, but  a completely different country. I don't actually have any evidence bats are eaten frequently in China.


text:
Interpretion of my coronavirus 2019-nCov, Wuhan, China BLAST tree?
This is the BLAST tree of the latest coronavirus out of China (from Wuhan Institute of Virology, China). It seems strange that there is so much divergence from all the other coronaviruses. Is this expected of new diseases?

the key is for me how much divergent and where (someone speak about that on Twitter) if the divergence is only in a spike protein and in some region, is that natural possible? what can explain all? I heard one bat virus has protein E match 100% with sars2 envelope protein... the only thing can explain all this is one ogm virus...

The  evolutionary related group (clade) of betacoronaviruses you have identified share an amino acid homology of 85% and include SARS. I know this from the underlying tree published on BioRxiv of a broader group of betacoronaviruses, i.e. your data is a defined subset of the betacoronaviruses which all share a unique, single common ancester.
Lets call this group the SARS clade.
You have performed a nucleotide blast and asked NCBI to produce a NJ tree using 2019-nCov as the reference. I can tell this by the colour coding and the scale bar at the bottom right hand corner shows the genetic distance involved is reasonably larger than 15% divergence by amino acid data. The scale bar represents the number of mutations per nucleotide.
In summary, your tree is essentially a subset of the broader genetic diversity of the one published, but there is a rooting issue.
In your tree, the majority of sequences are from the 2002 SARS epidemic and the virtually zero genetic distance between them is simply because it is a rapidly transmitted outbreak. I didn't realize that SARS has two independent origins both initially from bats. This is quite scary.

2019-nCov is an outgroup within the 'SARS clade' hence it appears at the other side of the tree, i.e. they share a more distant common ancestor
However 2019-nCov is not the most distant common ancester within this subset of the betacoronaviruses, i.e. SARS clade, this belongs to the two bat viruses ZC45 and ZX21. The program has made a likely rooting error (below).
Again the blast omitted the majority of the betacoronaviruses, for example MERS

Rooting issue the reason I suggest there is a rooting error is because the BioRxiv tree using a broad sample of the betacoronaviruses placed the bat strains ZC45 and ZX21 as outgroups to the SARS clade and 2019-nCOV was immediately within that, so in this definition 2019-nCOV is an 'ingroup' within the SARS clade, whereas your tree it is an 'outgroup'. It not a huge issue, but the location of the "root" is determined by the common ancester above it (sequences with <85% homology) and in your tree those have been omitted.
Generally, I like the analysis otherwise and gives an insight into SARS that I hadn't previosly been aware of.


text:
Clustered column chart for single cell data
I want to create a column plot for my single cell analysis that has 2 sample types, a normal and PDAC. I want to show, the percentage of different cell types per sample type.
How do I go about determining the percent of each cell type in each condition and then go about plotting the column plot? Also, is it possible to determine the SD from the dataset and add it?
dput(head(all@meta.data))
structure(list(orig.ident = c("PDAC", "PDAC", "PDAC", "PDAC",
"PDAC", "PDAC"), nCount_RNA = c(7945, 7616, 7849, 5499, 853,
1039), nFeature_RNA = c(2497L, 2272L, 2303L, 2229L, 509L, 588L
), percent.mt = c(2.63058527375708, 3.7421218487395, 4.9433048796025,
0.490998363338789, 0.234466588511137, 0.192492781520693), RNA_snn_res.0.5 = structure(c(5L,
5L, 5L, 5L, 1L, 6L), .Label = c("0", "1", "2", "3", "4", "5",
"6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16",
"17", "18", "19"), class = "factor"), seurat_clusters = structure(c(5L,
5L, 5L, 5L, 1L, 6L), .Label = c("0", "1", "2", "3", "4", "5",
"6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16",
"17", "18", "19"), class = "factor"), S.Score = c(0.0171539298241524,
-0.00599651468836135, 0.011058574403656, -0.0424413740104151,
0.03549521163095, -0.0243347616753686), G2M.Score = c(-0.0591387992008088,
0.0253275642795205, -0.0512402869839816, -0.0239076248512967,
-0.0263164126515597, -0.0339086517371782), Phase = structure(c(3L,
2L, 3L, 1L, 3L, 1L), .Label = c("G1", "G2M", "S"), class = "factor"),
    old.ident = structure(c(1L, 1L, 1L, 1L, 1L, 3L), .Label = c("Fibroblast",
    "T Cell", "Endothelial", "Tumor", "Stellate", "Macrophage",
    "B Cell", "Mast Cell", "Acinar", "Endocrine", "Exocrine"), class = "factor")), row.names = c("PDAC_AAACCCAGTCGGCTAC-1",
"PDAC_AAACGAAGTCCAGGTC-1", "PDAC_AAAGAACAGCAAGTGC-1", "PDAC_AAAGAACAGGATGAGA-1",
"PDAC_AAAGTGACATCCGTTC-1", "PDAC_AACAACCAGGAAGTAG-1"), class = "data.frame")

You can start with prop.table() to get the percentages.

@haci So I would do prop.table(all@meta.data)

prop.table() requires a table as input so you can try something like prop.table(table(object@meta.data$column1, object@meta.data$column2).

You can make use of table() and prop.table() to calculate the number and proportions of your cell types with respect to tissue of origin.
prop.table() requires a table as input so you can start with something like prop.table(table(object@meta.data$cell_type, object@meta.data$tissue). This will create an objet of class table.
Next step is to create a data frame with as.data.frame(your_table_object) and then to convert it to "long format", for example with tidyr::pivot_longer(). You can feed these data in long format to ggplot().

If I understand your question correctly and you have dplyr and ggplot installed, you could try something like (if a is your data.frame)
a %>%
dplyr::group_by(orig.ident, old.ident) %>%
dplyr::summarise(count  = n()) %>%
ggplot(aes(x=orig.ident, y=count)) + geom_bar(stat='identity', aes(group=old.ident), position='dodge')


text:
How to represent pathways as mathematical graphs?
I would like to construct a mathematical graph: $G = (E,V, w(E))$ a graph with vertices, edges with weight on them based on a KEGG pathway.
Is there a way to do that with existing tools? If not I will start to implement a Python module for that. My main question is what is the mathematical graph representation of pathways? Assuming the possible node types and the different connections?
I would like to transform a pathway to a graph because then one can run min-cut-max-flow algorithms and stuff like that as part of my pathway analysis.

Pathways are generally directed graphs, which can be implemented in python using the networkx package. Note that there will be cycles in graphs.
Please note that there are large numbers of pathways analysis packages that already exist (both python and R packages but also just webservers), so I strongly recommend searching to see if something already exists before wasting time and implementing it yourself.

Pathways aren’t exactly directed graphs. There are different relationships (inhibition, enhancement,...) how to build a netwrokx’s graph out of it isn’t clear. And that is my main issue, how to build a graph given a KEGG pathway.

Basically my original question is how to transfer a KEGG pathway to a networkx object?

Inhibition and enhancement would be weights with different signs. For more complicated interactions standard graph representations simply wouldn't suffice.


text:
Kallisto for Nanopore?
Although the Kallisto software was designed for RNA, it has also been very reliable in counting selected FASTA sequences in DNA sequencing. For example, given a  WGS sample and my FASTA file. Kallisto is able to accurately return the counts for each sequence in my FASTA file.
Now, I'd like to do something similar on Nanopore long reads. Is that a similar software?

It seems NanoCount can do what you are looking for.:

EM based transcript abundance from nanopore reads mapped to a
  transcriptome with minimap2 Python package adapted from
  https://github.com/jts/nanopore-rna-analysis by Jared Simpson
NanoCount estimates transcript abundance from ONT direct-RNA
  Sequencing reads mapped to a transcriptome. It uses an
  expectation-maximization approach like RSEM, Kallisto, salmon, etc to
  handle multi-mapping reads. The reads must be mapped to the transcript
  set using minimap2.

Thanks. But the page (https://github.com/a-slide/NanoCount) doesn't show an option for me to provide a FASTA file. How does the program know what transcripts I have in my example?

For example, ENSDART00000010144 transcript name. How did that go inside the program (as reported on the page)? Sorry I just don't get how it'd work without asking the user to index like Kallisto.

According to the documentation it expect that you first align your reads to a transcriptome with minimap2, and provide that BAM or SAM file to NanoCount.


text:
How to obtain file offsets corresponding to range query in block-gzipped file
How can we retrieve the file offsets of the zipped-chunks in a block-gzipped file that contain records for a region/window of interest (e.g. Chr01:100-Chr02:145)?
It seems like the CLI of samtools, or tabix should support that, but I can't find a way to just get the offsets in the datafile, without having the accompanying datafile at hand. I couldn't find a formal specification for those index files either.
Perhaps one can recommend a library that can parse those index files, in python, C, or golang? It's hard to accept that there isn't a reusable block-gzipped index parser somewhere, other than the one buried deep in samtools/htslib/gatk code.
edit
I'd initially given a small motivating example making use of file offsets to access a subportion of file on a remote server. Commenters were quick to point out that samtools supported a series of protocols for doing exactly that: ftp, http, and s3 (which requires https, which requires that samtools be compiled with libcurl support) Thanks! That example was detracting, so I removed it.
samtools view does fetch a subportion of a file (potentially from a remote location), but it does a few extra things, like decompressing the chunks retrieved. With this question, I'm interested in manipulating the gz chunks directly.

samtools and tabix support remote ranged file queries, so why bother implementing it yourself?

It supports remote ranged queries only for ftp: http: I picked a poor use case example. I actually need to access the file over a protocol not supported by samtools -- S3.  It is https, but has custom authentication headers. I might also want to do the same over ssh or some other protocol, or via custom authentication scheme, like ntlm or whatnot. I was trying to avoid making this an s3 or protocol question -- I just need the file offsets.

Here is the full documentation: https://www.htslib.org/doc/htslib-s3-plugin.html

ugh. they're going to support every protocol under the sun eventually? There's already very good s3 clients out there. Will it support datasync, parallel uploads, separation of each gzip block into its own s3 object multipart fragment too to take advantage of better replication? Strange decision imo. What if my bam file's split in a million files?

It's supported via a plugin, so presumably there will just be more of those in the future.

So if the purpose of obtaining file offsets is to access subportions of remote files via protocols that samtools supports, as mentioned in this comment, then you can make your life easier and use samtools/tabix directly on the cli, prefixing your datafile with the corresponding scheme.
One side effect of using the CLI however, is that the blocks retrieved will be decompressed (e.g. samtools view). You'd have to re-compress them if what you're trying to do is create a smaller "bam" file.
Keep in mind that if you retrieve all zipped chunks which cover your range of interest, then you might be extracting more data than you asked for. The first and last chunk may contain records outside your region, just because their neighboring records happened to be compressed in the same chunk. You'd have to be careful to filter those downstream.
Specification for the tabix format
If you absolutely need a custom solution, the tabix data structure format, authored by Heng Li, is on the samtools website:
https://samtools.github.io/hts-specs/tabix.pdf
A more in depth discussion of how the data structure works can be found here:
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3042176/
The tree-shaped binning index is a bit hard to understand from the text alone, but the tabix paper mentions that their approach is almost identical to the binning algorithm in this paper: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC186604/pdf/X6.pdf (page 1003). That paper has a nice illustration which makes things clearer.
Operation
Note: This is from my understanding of how the indices interoperate.
You don't need the data file to retrieve chunk offsets, you can do the work with just the tabix file. But it's complicated, and probably error prone if you're implementing this yourself.
Let's assume the query interval is Chr01:100-Chr02:2000.
Each sequence (Chr01, Chr02, etc) has its own pair of indices (binning index + linear index) in the tabix file. So you would break the query into two queries: Chr01:100-999999999999 and Chr02:0-2000, and concatenate the two results.
Within one sequence, to find all the blocks for open interval [X,Y), you would use the binning index to retrieve all the bins overlapping this interval (There is C code to do this in the data structure pdf). Then, you need to go through all the possible bins selected, and filter chunks. Not all chunks in a given bin will be overlapping with X, Y.
To filter chunks within a bin, you need the help of the linear index. There is one entry in the linear index per 16kbp of sequence positions. Each entry of the linear index contains the virtual offset of the leftmost record overlapping that 16kb window. I don't know what the entry contains if there is no record overlapping that 16kb window, presumably 0.
You would find the linear entry in which position X falls (X/16Kb), take that offset, call it O_x. You find the linear entry corresponding to position Y, and take that offset, call it O_y. When you go through the selected bins, filter out any chunk whose offset is < O_x, or whose offset >= O_y.
Assuming your datafile is sorted, the chunk range you're interested in is contained between the minimum offset and the maximum offset selected.
The virtual offsets are 64bit unsigned integers. The higher 48 bits of the virtual chunk offset corresponds to the file offset of the datafile (the other 16 you don't need. They are the offset of the record within that decompressed chunk).
Software
The software I have found typically mimics closely the hts_lib API. And hts_lib doesn't export a function which loads an index without also requiring a datafile. As far as I know.

From the htslib side of things there'd be no reason to load an index without a corresponding data file, since there are already functions for performing the range queries.

For htslib, the APIs you would be interested in are those discussed in htslib #385. For various reasons, something similar will likely make it into HTSlib in the not too distant future.


text:
Converting values to string
I have a data frame
> head(custom.cn.data)
      Gene       Sample_name         CN
1: DDX11L1          F02         0.6788844
2: DDX11L1          C02        -0.9924524
3: DDX11L1          A07        -2.3833559
>

I want to convert values in CN column < 0 to Del and >0 to Amp
But I am getting error
> a %>% mutate(q = replace(CN, which(CN<0), Del))
Error in replace(CN, which(CN < 0), Del) : object 'Del' not found

> a$CN <- replace(a$CN, which(a$CN < 0), Del)
Error in replace(a$CN, which(a$CN < 0), Del) : object 'Del' not found

Any help?

You are getting these errors because you did not use quotes before and after Del and R looks for an object called Del instead of using the string "Del".

What happens to when CN is zero?

custom.cn.data$new_column <- ifelse(custom.cn.data$CN < 0, "Del", "Amp")

In this scenario CN == 0 would be Amp

The following code would create a new column with the desired values, 0 values will have NA in this new column.
custom.cn.data[custom.cn.data$CN < 0, "CN_string"] <- "Del"
custom.cn.data[custom.cn.data$CN > 0, "CN_string"] <- "Amp"

Try to stay within one package (paradigm) for one problem. This will make the code more readable, easier to debug.
At the moment you have data as a data.table object from data.table package, then trying to add a new column using mutate from dplyr package, and within mutate using base replace function.
Solutions:
data.table
# update by row index
custom.cn.data[ CN < 0, NewColumn := "Del" ]
custom.cn.data[ CN >= 0, NewColumn := "Amp" ]

or use ifelse
custom.cn.data[, NewColumn := ifelse(CN < 0, "Del", "Amp") ]

base
see @DavyCat's and @haci's answers.
dplyr
Note, this will change the data.table into data.frame.
# ifelse
custom.cn.data <- custom.cn.data %>%
  mutate(NewColumn = ifelse(CN < 0, "Del", "Amp"))

or use case_when
custom.cn.data <- custom.cn.data %>%
  mutate(NewColumn = case_when(
    CN < 0 ~ "Del",
    CN >= 0 ~ "Amp"))


text:
How to calculate coverage depth of genes in an annotated sequence?
I have a metagenomic dataset across multiple samples that has been grouped and binned (we did this due to poor sequence quallity).
It is several samples taken from termite guts of differing castes, the samples have been grouped to make binning easier due to poor sequence suitability. They have been binned via ASVs, assembling reads from all samples to generate the putative bacterial species. So I have fasta files of putative species that have been annotated via prodigal. I want to now calculate the abundance of each bin in a given sample. I think the best way to do this would be at the gene level by looking at the number of reads that map to a gene in any given sample. I am just unsure how to do this.
I have annotated the bins and now want to calculate the abundance of each bin in a sample but looking at the number of reads from a sample that map to each gene in a bin.
Is this a viable method? If so how would I do it and if not what else could I do?

The best way to approach it would be to take trimmed reads map them to the bins (i.e. with bowtie2) and estimate the abundance of each bin with samtools by getting the depth of coverage:
https://www.metagenomics.wiki/tools/samtools/breadth-of-coverage
or you can use QUANT_BIN module from the METAWRAP to do it for you. Follow this for the METAWRAP approach:
https://github.com/bxlab/metaWRAP/issues/40


text:
samtools / bamUtil | Meaning of <wildcard> as Reference Name
I have been working with several tools on bam files recently and I am not sure how I should interpret some of the outputs.
samtools idxstats mapped.sorted.bam

Chr1    15096   0   0
Chr2    33397   0   0
Chr3    43888   41  0
Chr4    20819   1   0
*   0   0   34

Documentation says third and fourth column give "mapped reads" and "unmapped reads". So I interpreted the last line (*) as the number of unmapped reads as this is the only line with a non 0 value at this field.
Then I used
bam splitChromosome --in mapped.sorted.bam --out test/aln.

Reference Name: chr6 has 6 records
Reference Name: Chr7 has 1 records
Reference Name: Chr3 has 41 records
Reference Name: Chr4 has 1 records
Reference Name: * has 34 records

And looking at the bams outputted (one per reference sequence) in test/ I find an aln.UnknownChrom.bam that I am not sure where it comes from. I actually found 97 unique reads ID in it, where I was waiting for 34 or less.
So my first question is : Does "*" actually means "all reads that didn't mapp" and this information is stored in way I didn't know about in BAM files ?
And my second question would be : Anybody knows where this aln.UnknownChrom.bam comes from when using splitChromosome ? It is a source of error in a pipeline I am building and I'd like to get rid of it, if these are not actual mapping informations.

As per the SAM/BAM format specifications:

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

So, yes, these are the unmapped reads.
I haven't used this splitChromosome utility before myself, but as it splits the BAM file into one BAM file per chromosome, I suspect it simply also creates a BAM file for the unknown (*) chromosome. Why the numbers don't match I can't tell.
How to deal with that file in your pipeline will depend on what the pipeline is meant to do. But if you aren't interested in these reads, I suppose you could simply remove it before continuing to the next step.

Nice ! I'll maybe ask directly on bamUtil's github for my second question.


text:
Filtering out 8-oxoG-derived artefacts in variant calls from cfDNA
I have cfDNA samples which were sequenced (at about 300-1000X) on a 30-gene panel via a capture-seq protocol for the purposes of calling variants on circulating tumour DNA.
However, the resulting variant calls (from both varscan and Mutect2) are dominated by transversions from C>A and G>T at low allele frequencies, which are characteristic of DNA oxidation of guanine to 8-oxoguanine, and which I suspect are artefacts.
I have matched baseline tumour samples from the same patients, sequenced in the same way, (and there are a few C>A/G>T transversions at low AFs in these samples that I can filter out by a direct allele fraction filter).
My question is this, what is the best way of filtering out the false positive mutations derived from the 8-oxoG formation and retaining the true positives (which will be at low allelic fraction)?
Some background: I do have serial (cfDNA) time points for some patients, so I am loath to just look at mutations called in the baseline tumour samples, as any (true) novel variants appearing later in the time points would be missed.


text:
Make a fasta out of mapped reads without taking into account the reference
I've aligned reads onto an haploid reference genome. I'd like to have a consensus sequence of all my aligned reads that doesn't take into account the reference genome I used (i.e. I just want the nucleotides from my reads at each base (taking into account frequency and quality) but if there is no read a one position I rather have N/? than the reference nucleotide.
Is there a way to do this ? Thanks

vSNP outputs a VCF file that includes positions with no coverage, represented as "N".  Then there is a tool to merge a VCF file into the FASTA reference  vsnp_merge_vcf_into_fasta.py.  It will merge those "N" regions into the FASTA reference. There are options for considering frequency and quality.
 It's limited to just SNPs and those regions along the reference with no coverage.  It does not merge small indels or unmapped reads.  Pilon is a better tool for applying indels, but I'm not aware of it letting regions with no coverage be represented as "N".
vSNP_step1.py -r1 *_R1*fastq.gz -r2 *_R2*fastq.gz -r ref.fasta

Use *zc.vcf output to merge into FASTA
vsnp_merge_vcf_into_fasta.py -f ref.fasta -v *_zc.vcf

This is similar to this question, which links to this solution.
The difference here is that you want to go through and replace low coverage (maybe zero-coverage?) positions with Ns rather than using the reference. You can likely accomplish this using samtools depth to measure coverage and then writing a script that will go through each position in your reference, check its coverage, and string-replace with an N. Probably a bash one-liner could do this if someone very devoted to bash one-liners came up with one (I don't care for them, but other people are more creative than I). Very quick and dirty pseudocode might look like this:
### Someone could come up with something better if they spent >3min on it.
### should be easy to get it working if you use biopython
# read_and_preprocess() stands for whatever you need to do to get data in the right shape
coverage_data = read_and_preprocess('samtools_coverage.txt')
# reference_with_variants.fasta comes from the solution linked above
reference_fasta = read_and_preprocess('reference_with_variants.fasta')

threshold = 0
fasta_string = ''
for position in reference_fasta:
  if coverage_data[position] <= threshold:
    fasta_string += "N"
  else:
    fasta_string += reference_fasta[position]

# nonexistent fn that writes a fasta.
write_as_fasta(fasta_string)

I am a little skeptical of this approach because it seems like if you do not trust the reference sequence to stand in for zero-coverage areas, then you probably shouldn't trust it to align to and call variants against. Read mapping can be a rather noisy process, especially when you are using a more distantly related reference. I think this is related to the @swbarnes2 comment.
I would suggest instead indicating the regions of low/zero coverage by putting them in lowercase acgt with your covered regions uppercase ACGT, which preserves the coverage information that you care about but still allows you access to the reference information. It's not a perfect solution, but I think it makes the most out of the marginal data you're working with, and it is trivial to change the above example code to do it.
Hope that helps-

Have you googled de novo assembly?

I can’t do a de-novo assembly as I have really low coverage. I’m working with really short reads from ancient DNA.

You can make a consensus fasta from your bam, but that uses the reference.  If your read depth is that low, can you be confident of discrepancies between your reads and the reference?

No, but that’s another problem. Thanks. The thing I want to do can also be useful when you want to reassemble only some targeted genes from medium coverage data I guess.


text:
Is it possible to do this in R?
I have a big clinical data file with 96 columns like age, gender, BMI, etc
I want to see which of these clinical characteristics respond to chemotherapy. Response to chemotherapy divides patients to two groups as yes and no. Among these clinical characteristics, some are based on continuous data and some reflect `categorical data.  Is it possible to write a function to loop over these characteristics and test for relationship between response to chemotherapy? Something like Wilcox test for continuous values and chi square for categorical characteristics. By hand, it takes too much to do for all one by one !
This is my data
ID  response_to_chemo   BMI Chemo   Predictor   DJANGO  gender
AH/155  no  50  1   1   0   M
RS/022  no  67  1   1   0   M
RS/027  no  80  1   1   1   F
ST/023  no  65  1   1   0   M
SH/051  yes 47  1   1   0   M
AH/075  yes 90  0   1   0   M
RS/047  no  67  0   1   0   F
ST/029  no  61  1   1   0   F

For instance
myTable <- table(clin$gender, clin$response_to_chemo)

chisq.test(myTable)

Gives p-value for testing the relationship of gender (CATEGORICAL) in response to chemotherapy
And
t.test(clin$BMI ~ clin$response_to_chemo)

is for BMI
I meant a function calculating these p-values for clinical characteristics one by one

Did you triy anything yet ?

Actually in my mind I think about some for i in colnames of my data but can not go further as I am unable to tailor a function :(

Try to include some code that you would run on the data example you provided. Taking the time to make a good reproducible example often gives you the answer to the original question.

Don't forget to correct for multiple hypothesis testing at the end, too. With 20 completely random variables, you'd expect 1 of them to show a spurious association at p<0.05.

R supports logistic regression, which would seem to be the most efficient method for tackling this question. Assuming the "Chemo" variable is the type of chemo the  code would be something like:
glm( (response_to_chemo == "yes") ~  BMI + Chemo  + Predictor +  DJANGO + gender,
        family="binomial")

EDIT: Corrected a typo (added a missing double quote)

The general approach would be to loop through every column starting at column 2. You can use numeric indexes to do that.
For each column, check its type. If it is a factor, use your chisq.test method. If it is of numeric type, use your t.test approach. Write the results to a list so you can parse them later.

You have category-based variables as well as variables based on continuous data. The best-option could be linear Mixed- effects model for your research purpose. Chi-square test and t-test may be used to cross-validate or re-examine the outcome of LME  model.

Hi @subhash what mixed model would you consider?

When using such a model, how do you decide which effects to be "fixed" and which to be "random"?

These are assumptions that pervade hierarchical modelling.Saying "which effects" here seems to reflect the variables! Please clarify.


text:
How to get the splicing regions (BED) of all genes for Human GRCh38?
How can I obtain a BED file that describes the splicing sites of all genes for GRCh38 human genome ?
Thank you

Interested to know if there are special file format for RNAsplicing.

Is it the alternative splicing events you are looking for? https://bioconductor.org/packages/devel/data/annotation/manuals/alternativeSplicingEvents.hg38/man/alternativeSplicingEvents.hg38.pdf

@PaulEndymion thanks for the link. Does the alternative splicing events mean all splicing regions, even for the same gene ? if that's the case, then yes, I am looking for something like that !


text:
How to exclude the repetition of gene-gene correlation calculation in python?
I am trying to make a paired matrix of gene-gene correlation. Considering that I have a huge matrix (13000 genes and 900 samples) and for some reasons I don't want to decrease the number of my genes, my gene-correlation matrix would be 13000*13000 and my paired matrix will become 169 million *4 (Column 1: Gene 1; Column 2 : Gene 2; Column 3: Correlations; Column 4: P-values) . In this case, I have to exclude unnecessary calculations as much as I can. I have excluded the situation that Gene 1 = Gene 2. But I couldn't find a way to exclude the condition that "Column 1: Gene 1 ; Column 2: Gene 2 = Column 1: Gene 2; Column 2: Gene 1 ". To make a long story short, correlation between G1 and G2 is equal to G2 and G1. It is like calculating just lower section of diag in a symmetric matrix. I would be grateful if anybody help me in this case. I have enclosed my python codes here:
import pandas as pd
import numpy as np
import scipy
import math
import openpyxl
from openpyxl import Workbook
from scipy.stats import spearmanr

din=pd.read_csv('m_test.csv', index_col=0)

out=pd.DataFrame()
outdf=pd.DataFrame()

for g1 in din.index:
    for g2 in din.index:
        temp=din.loc[[g1, g2]]
        if g1==g2:
            next
        else:
            spR, spP=spearmanr(temp.loc[g1], temp.loc[g2])
            frame={'g1':[g1], 'g2':[g2], 'spR':[spR], 'spP':[spP]}

            out=pd.DataFrame(frame)

            outdf=pd.concat([outdf, out])

Actually I found an easy and efficient way.
import pandas as pd
import numpy as np
import scipy
import math
import openpyxl
from openpyxl import Workbook
from scipy.stats import pearsonr
from scipy.stats import spearmanr
from scipy.stats import kendalltau

pro=pd.read_csv('m_test.csv', index_col=0)

xInds = []
yInds = []
for i in range(len(pro.index)):
    for j in range(i+1, len(pro.index)):
        xInds.append(i)
        yInds.append(j)

z=0
Rlist_sp = []
Plist_sp = []
Rlist_pe = []
Plist_pe = []
while z < len(xInds):
    b=xInds[z]
    c=yInds[z]

    spR, spP = spearmanr(pro.iloc[b].values, pro.iloc[c].values)
    peR, peP = pearsonr(pro.iloc[b].values, pro.iloc[c].values)

    Rlist_sp.append(spR)
    Plist_sp.append(spP)
    Rlist_pe.append(peR)
    Plist_pe.append(peP)
    z=z+1

G1=pd.DataFrame(xInds)
G2=pd.DataFrame(yInds)
R_sp=pd.DataFrame(Rlist_sp)
P_sp=pd.DataFrame(Plist_sp)
R_pe=pd.DataFrame(Rlist_pe)
P_pe=pd.DataFrame(Plist_pe)

Final=pd.concat([G1, G2, R_sp, P_sp, R_pe, P_pe], axis=1)
np.savetxt('Final_gg_A549.txt', Final.values, fmt='%s', delimiter='\t')

You don't need g2 to go from 0 to din.index, just from 0 to g1 - 1. This way, you'll end up calculating just for the "lower triangle":
g1-> 0   1   2   3
g2
0                    #(0,0 to 0,1-1 = nothing)
1    X               #(1,0 to 1,1-1 = 1,0)
2    X   X           #(2,0 to 2,2-1 = 2,0; 2,1)
3    X   X   X       #(3,0 to 3,3-1 = 3,0; 3,1 and 3,2)
4    X   X   X   X   #(4,0 to 4,4-1 = 4,0; 4,1; 4,2 and 4,3)

Thank you very much. It works this way.


text:
A new paper suggests the Corona Virus has "Uncanny similarity of unique inserts in the 2019-nCoV spike protein to HIV-1" - What does this mean?
Quote:

We found 4 insertions in the spike glycoprotein (S) which are unique to the 2019-nCoV and are not present in other coronaviruses. Importantly, amino acid residues in all the 4 inserts have identity or similarity to those in the HIV-1 gp120 or HIV-1 Gag. Interestingly, despite the inserts being discontinuous on the primary amino acid sequence, 3D-modelling of the 2019-nCoV suggests that they converge to constitute the receptor binding site. The finding of 4 unique inserts in the 2019-nCoV, all of which have identity/similarity to amino acid residues in key structural proteins of HIV-1 is unlikely to be fortuitous in nature.

Source: https://www.biorxiv.org/content/10.1101/2020.01.30.927871v1

I know this is a preprint paper and it's not yet peer reviewed, but can someone tell me what the implications are if this is true? Does this mean that the virus is artificially constructed?

UPDATE: The article has now been withdrawn with the following note:

This paper has been withdrawn by its authors. They intend to revise it
  in response to comments received from the research community on their
  technical approach and their interpretation of the results. If you
  have any questions, please contact the corresponding author.

This is very odd, and will require a rigorous investigation, but my initial reaction is one of scepticism. Considering just the 1st insert, the insert sequence is GTNGTKR, short at just 7 amino acids. A simple BLASTP vs NR did not find perfect matches to HIV sequences, but did reveal 100% identity across the full sequence length to >50 other short protein sequences, which could be spurious 'chance' hits of course. Many of the top 100 hits are against eukaryotic sequences, for example one is against Pristionchus pacificus, a type of nematode worm (see alignment below). Someone needs to do a proper peer review of this preprint before any conclusions are drawn from it.
>tank-1 [Pristionchus pacificus]
Sequence ID: PDM74036.1 Length: 2481
Range 1: 1474 to 1480

Score:24.0 bits(49), Expect:2477,
Method:,
Identities:7/7(100%), Positives:7/7(100%), Gaps:0/7(0%)

Query  1     GTNGTKR  7
             GTNGTKR
Sbjct  1474  GTNGTKR  1480

UPDATE 1: Things are moving fast. There are now 10 comments under the preprint, all of them critical of the idea that these inserts are meaningfully similar to HIV-1 sequences for the reason I outlined above (i.e. short sequences with many hits against organisms from across the tree of life) and additional points, such as small insertions being quite normal evolution for RNA viruses. There are also similar critiques on twitter and biorxiv has added the following disclaimer header to their website:

bioRxiv is receiving many new papers on coronavirus 2019-nCoV.   A
  reminder: these are preliminary reports that have not been
  peer-reviewed. They should not be regarded as conclusive, guide
  clinical practice/health-related behavior, or be reported in news
  media as established information.

I have never seen such rapid (and unanimous) post-publication peer review before. It is therefore already (just 1 day after publication) fairly clear that at least some claims in this study are completely false and there is no good evidence to support the claim that the 2019-nCoV have acquired sequences from HIV or evidence to suggest the virus is engineered.

Thank you for this. It really makes me wonder why a respected scientist would publish such an obviously flawed paper though?

Blast isn't the issue here, two proteins can been unrelated by blast but retain the same function. Functionally these proteins have similarities. The universal role of indels or what the authors call "inserts" particularly in viral surface antigens is widespread and has never been resolved. It is wrong to say the paper is "flawed", I think the observation is spurious, but has value.

@SurpriseDog I don't know the motivations of the authors but the article has now been withdrawn

Normally "inserts" used in the manuscript are "indels" in protein alignments, short for insertions and deletions.
What I think has happened is a group investigating indels in HIV env noticed indels in 2019-nCov. Essentially I think the correlation is spurious - but I haven't test it, but the area of research in understanding indels is certainly valid and important.
What is certain is that indels induce a large structural change to a protein structure and any Gibbs free-energy style calculation will identify this.
Vaccine
The spike protein will be the primary candidate to make a 2019-nCov vaccine and that is a very important reason why the sequence was rapidly released. So it is an important protein and the structural changes indels induce mean that a SARS vaccine will probably not provide much protection against 2019-Cov, even apart from the amino acid divergence (below).
Differences 2019-nCov vs HIV
In summary, alot. HIV env and particularly HIV gag are very different from coronaviruses, both in the mechanism of genome replication, coronavirus never leaves the cytoplasm, clinical outcomes, tissue tropism and duration of infection.
Similarities
HIV env and the glycoprotein spike of coronaviruses are the receptor binding protein to gain entry into a cell. They are called structural proteins. Entry to a cell can be blocked by antibodies and these antibodies are called "neutralizing antibodies". Neutralizing antibodies are catastrophic for a virus. Other antibody responses can be effective, such as IgM, but to clear an infection just using antibodies, you need neutralising antibodies. Both HIV env and the coronavirus spike protein are subject to neutralising antibodies. HIV gag has nothing to do with HIV env, in terms of function or antibody exposure. This is why the spike protein will be the primary vaccine candidate for a subunit vaccine.
Coincidence, law of chance
There is large variation of indels in HIV env within HIV and what the authors are inferring is there is a resemblence to that between SARS and 2019-nCov. In my opinion this is a coincidence, because they are comparing a large repertoire of HIV varients against a single indel pattern in the coronaviruses.
Why coronavirus indels?
That is a very good question. Generically indels in viral surface antigen genes are common, much more common in other proteins - such as those involved in virus replication (non-structural proteins). The amino acid identity between SARS and 2019-nCov is 80%, and in any virus, such as flaviviruses 80% identity means indels will be present in surface antigens between the viruses. The answer is it is not unusual in any RNA virus to see indels at a comparatively large amino acid divergence.
What function could they serve
I've briefly looked at indel bioinformatics between flaviviruses (Zika virus, yellow fever virus etc..) notably using envelope (E) protein sequences, and they also occur between African Zika viruses in the E-protein. E-protein being the equivalent of coronavirus spike protein, the receptor-binding protein. No-one has ascribed a function to them and that is the problem with this manuscript.
Hypthoses

One theory is that a structural change in the protein will occur to stop antibody binding.
Another theory is they have functional differences, such as cell tropism

Bioinformatically separating the two theories is extremely hard without wetlab experimentation.

Really thorough explanation, clear and precise !! I learned a lot thanks ;)

"Neutralizing antibodies are catastrophic for a virus. " why because target the T cells to the virus? and if the virus can enter in T cells destroy them and reproduce too what happen???

Good common. No  HIV targets T-cells. I will update the post when I have time "neutralising" means preventing the virus binding to the cellular receptor, thus it cannot enter a cell and therefore cannot replicate

There's a new paper by a Chinese Virologist who worked at a WHO reference lab and has now fled to America claiming whistle blower status. Here is the paper: https://zenodo.org/record/4028830#.X19xByXZglR - Any thoughts? I'm hesitant to post another question about stuff like this, being so far out of my depth.

I'll have a look when I've a moment

SurpriseDog - I suspect the wording "unlikely to be fortuitous in nature" is what led you to ask "Does this mean the virus is artificially constructed?" - No!
I think the writing here is imprecise and confusing. There is no implication that these are human-generated (non-natural/engineered) insertions but rather I think they are implying that these mutations are likely to have functional implications. Whether there is any relation to HIV is questionable.
There has been some political finger pointing and threats of full-scale investigation into whether the virus escaped from a lab. Hard to say, but current thinking is that SARS-Cov-2 is a natural isolate, not engineered (see https://www.nature.com/articles/s41591-020-0820-9). It was probably transmitted from bats.


text:
How to generate conformations for a ligand for in silico screening?
I have multiple ligands (in pdb format) which I want to try and dock to various proteins.
But for each ligand I have only a single conformation (by which I meant the spatial relationship between it's atoms, not its position relative to the target protein).
How do I generate other possible conformations?
Is there a generally used tool for that?

Still not clear what you mean by “pose”. You need to use standard scientific English on this site.  Do you mean structure? If not it shouldn’t be difficult to find a scientific paper that uses the correct term.

Thanks.
I see that these are properly referred to as conformations.

Ok. I assume that your ligands are peptides or proteins that can adopt different conformations. There may be software for doing this, but it will be quite sophisticated. Good luck.

I'm going to migrate this because it's a better fit for Bioinformatics since it's about tools rather than the underlying biology, and I think you're more likely to get a good answer there. I'm only loosely aware of their policies there so you may need to make edits.

@David "Pose" is standard scientific English for macromolecular docking.

@canadianer — I stand corrected, although wonder how many on this list are familiar with the term.

One option I found was Frog2
https://mobyle.rpbs.univ-paris-diderot.fr/cgi-bin/portal.py#forms::Frog2
https://github.com/tuffery/Frog2
From an input of a mol2 or sdf file it creates a plethora of conformers, which can be outputted as a pdb file

Open babel
The most commonly used tool is open-babel (obabel). It seems annoying that you have to specify the number of conformers, but this is because if you get something like retinol, the number of conformers is huge. The command from the documentation is:
obabel startingConformer.mol -O ga_conformers.sdf --conformer --nconf 30
   --score rmsd --writeconformers

An alternative argument is --confab which will pick each representative of a RMSD cutoff cluster.
Do note that 2D representations and un-protonated molecules will need extra command line arguments.
RDKit
You can also generate conformers in RDkit, but I think this is beyond your requirements, but I thought it may be good to mention the essential compchem python library.


text:
Bedtools wrongly indicates zero coverage
I have a BAM file produced using BWA-MEM and I need to calculate the coverage of each exon in this sequencing. I tried using this command:
bedtools coverage -a exome_sorted.bed -b exome_alignment.bam -sorted
And the result was the following:
chr1    12595   12802   ref|DDX11L1,ref|NR_046018,ens|ENST00000518655,ens|ENST00000450305,ens|ENST00000456328,ens|ENST00000515242  .       ref|DDX11L1     0       0       207     0.0000000
chr1    13163   13658   ref|DDX11L1,ref|NR_046018,ens|ENST00000518655,ens|ENST00000450305,ens|ENST00000456328,ens|ENST00000515242  .       ref|DDX11L1     0       0       495     0.0000000
chr1    14620   15015   ref|WASH7P,ref|NR_024540,ens|ENST00000538476,ens|ENST00000438504,ens|ENST00000488147,ens|ENST00000541675,ens|ENST00000423562       .       ref|WASH7P      0       0       395     0.0000000
...

That is, zero coverage for all exons. But there are reads mapping to the exons, as I can see by opening my BAM file on IGV and looking at some exonic regions:

I have already tried sorting the BED file using sort -k1,1 -k2,2n (as recommended on the Bedtools documentation) and using bedtools sort with the -faidx option, but neither option solved my issue.
if I skip the -sorted option, the memory usage becomes way too high.
I have also checked if the chromosome names are the same in the BAM file and in the BED file, using samtools view -H. Both use chr1,chr2,...,chr22,chrX,chrY, although BAM has chrM and stuff like chrUn_gl000214 which are not on the BED.
Thanks!
EDIT:
@ATpoint pointed out that all visible reads in this region have MAPQ0. I looked at another exon that did not have this issue, but Bedtools still outputs zero coverage. IGV picture below:

If this is your output :
chr1    12595   12802   ref|DDX11L1,ref|NR_046018,ens|ENST00000518655,ens|ENST00000450305,ens|ENST00000456328,ens|ENST00000515242  .       ref|DDX11L1     0       0       207     0.0000000

Then I guess your bedfile looks like this (the first 6 fields) :

chrom : chr1
FeatureStart : 12595
FeatureEnd : 12802
name : ref|DDX11L1,ref|NR_046018,ens|ENST00000518655,ens|ENST00000450305,ens|ENST00000456328,ens|ENST00000515242
score : .
strand : ref|DDX11L1

The issue I see here is that you are giving the wrong information to the strand field of your bedfile. It should be either +, either - or . if you don't want to consider the strand information.
Hope it helps !


text:
Looking for a way to download all pdb files that detail a protein-drug complex
I am looking for a way to download all pdb files that detail a protein-drug complex.
While I could download the entirety of the pdb and search manually (or maybe write a script to discard all proteins that do not have an HETATM entry that details a non-water, non-salt molecule), that sounds effort intensive and error prone.
It seems reasonable that this subset already exist, but I was not able to find it.
Does such a list/database exists?

Albeit only a small fraction, the SGC has Target Enabling Packages which are few dozen sets of crystal structures each with several hits, i.e. non-natural targets, but not necessarily with nanomolar affinity but generally obeying Lipinski's rules (drug-like). For some times of analysis this wider but shorter set many be a better set than a messier set with only one target each...

Also, if you do end up looking manually, in this blog post is a list of crystallisation ligands, like glycerol, to blacklist.

Above I mention a small but diverse dataset with many ligands each. An alternative is KLIFS, which has only kinases, which is curated, has a lot of members and many ligands, but is not at all diverse. As a result it is a clean dataset often used as a testbed for novel methodologies. However, it represents only one family.

Proving a negative is hard, but to my knowledge I do not know a database that has all the structure of drug bound protein. I can give some close matches.
First a few terms to be on the same page. A ligand is something that binds, which can be physiological, such as a cofactor or substrate, or not (especially man-made), such as a drug, i.e. a compound that has a physiological effect, or a hit, a small fragment from a fragment screen that will not bind with nanomolar affinity. Generally, ligands are small molecules, but may include oligopeptides and oligonucleotides. In the case of crystal structures, you also have crystallisation additives, which are not ligands, but are hetatoms in PDBs. Also, the definition of drug may vary. drugbank.ca has the following categories of drugs "approved", "experimental", "withdrawn" for pharmaceuticals, "illicit" for narcotics and "investigational" for probes and counts macromolecules as drugs.

RCSB PDB: drug mapping table is a simple searchable table with ligands that are clinical drugs only —experimental compounds (=probes) will be absent.
bindingdb.org contains many table, including a '3D' one. This table contains all the PDBs with ligands and their affinities. This contains crystallisation additives and natural ligands. Peter Curran's repo for the HotSpots algorithm (first written by Chris Radoux) has a blacklist, but you are bound to get stuff missed.
bindingmoad.org is a database with all the protein-ligand structures from PDB up until 2018, but it includes natural ligands, like cofactors and substrates as above.
https://klifs.vu-compmedchem.nl a well curated database of kinases, which has a lot of members and many ligands, but is not at all diverse. As a result it is a clean dataset often used as a testbed for novel methodologies.
I mentioned in the comments the SGC Target Enabling Packages as a small subset of proteins with many hits each. But I realise that there is not a handy PDB list and the targets are not drugs, but probes.


text:
How to add an increment to a repeated strings
This is my data
 dput(METADATA)
structure(list(Run = c("SRR2920466", "SRR2920477", "SRR2920478",
"SRR2920506", "SRR2920507", "SRR2920532", "SRR2920531", "SRR2920469",
"SRR2920470", "SRR2920481", "SRR2920482", "SRR2920500", "SRR2920501",
"SRR2920536", "SRR2920537", "SRR2920471", "SRR2920472", "SRR2920483",
"SRR2920484", "SRR2920505", "SRR2920538", "SRR2920539", "SRR2920475",
"SRR2920476", "SRR2920487", "SRR2920488", "SRR2920542", "SRR2920543",
"SRR2920555", "SRR2920559", "SRR2920570", "SRR2920573", "SRR2920578",
"SRR2920583", "SRR2920586", "SRR2920594", "SRR2920546", "SRR2920547",
"SRR2920548", "SRR2920549", "SRR2920550", "SRR2920551", "SRR2920553",
"SRR2920558", "SRR2920563", "SRR2920565", "SRR2920567", "SRR2920569",
"SRR2920572", "SRR2920575", "SRR2920577", "SRR2920580", "SRR2920582",
"SRR2920585", "SRR2920589", "SRR2920590", "SRR2920591", "SRR2920593",
"SRR2920554", "SRR2920552", "SRR2920556", "SRR2920557", "SRR2920560",
"SRR2920561", "SRR2920562", "SRR2920564", "SRR2920566", "SRR2920568",
"SRR2920571", "SRR2920574", "SRR2920576", "SRR2920581", "SRR2920584",
"SRR2920587", "SRR2920588", "SRR2920592", "SRR2920595", "SRR2920579"
), source_name = c("HSC", "HSC", "HSC", "HSC", "HSC", "HSC",
"HSC", "CMP", "CMP", "CMP", "CMP", "CMP", "CMP", "CMP", "CMP",
"GMP", "GMP", "GMP", "GMP", "GMP", "GMP", "GMP", "Mono", "Mono",
"Mono", "Mono", "Mono", "Mono", "LSC", "LSC", "LSC", "LSC", "LSC",
"LSC", "LSC", "LSC", "Blast", "Blast", "Blast", "Blast", "Blast",
"Blast", "Blast", "Blast", "Blast", "Blast", "Blast", "Blast",
"Blast", "Blast", "Blast", "Blast", "Blast", "Blast", "Blast",
"Blast", "Blast", "Blast", "Blast", "pHSC", "pHSC", "pHSC", "pHSC",
"pHSC", "pHSC", "pHSC", "pHSC", "pHSC", "pHSC", "pHSC", "pHSC",
"pHSC", "pHSC", "pHSC", "pHSC", "pHSC", "pHSC", "pHSC")), class = c("spec_tbl_df",
"tbl_df", "tbl", "data.frame"), row.names = c(NA, -78L), spec = structure(list(
    cols = list(Run = structure(list(), class = c("collector_character",
    "collector")), `Assay Type` = structure(list(), class = c("collector_character",
    "collector")), AvgSpotLen = structure(list(), class = c("collector_double",
    "collector")), BioProject = structure(list(), class = c("collector_character",
    "collector")), BioSample = structure(list(), class = c("collector_character",
    "collector")), Cell_type = structure(list(), class = c("collector_character",
    "collector")), `Center Name` = structure(list(), class = c("collector_character",
    "collector")), Consent = structure(list(), class = c("collector_character",
    "collector")), `DATASTORE filetype` = structure(list(), class = c("collector_character",
    "collector")), `DATASTORE provider` = structure(list(), class = c("collector_character",
    "collector")), `DATASTORE region` = structure(list(), class = c("collector_character",
    "collector")), donorid = structure(list(), class = c("collector_character",
    "collector")), Experiment = structure(list(), class = c("collector_character",
    "collector")), GEO_Accession = structure(list(), class = c("collector_character",
    "collector")), Instrument = structure(list(), class = c("collector_character",
    "collector")), LibraryLayout = structure(list(), class = c("collector_character",
    "collector")), LibrarySelection = structure(list(), class = c("collector_character",
    "collector")), LibrarySource = structure(list(), class = c("collector_character",
    "collector")), MBases = structure(list(), class = c("collector_double",
    "collector")), MBytes = structure(list(), class = c("collector_double",
    "collector")), Organism = structure(list(), class = c("collector_character",
    "collector")), Platform = structure(list(), class = c("collector_character",
    "collector")), ReleaseDate = structure(list(format = ""), class = c("collector_datetime",
    "collector")), sample_acc = structure(list(), class = c("collector_character",
    "collector")), `Sample Name` = structure(list(), class = c("collector_character",
    "collector")), sample_type = structure(list(), class = c("collector_character",
    "collector")), source_name = structure(list(), class = c("collector_character",
    "collector")), `SRA Study` = structure(list(), class = c("collector_character",
    "collector"))), default = structure(list(), class = c("collector_guess",
    "collector")), skip = 1), class = "col_spec"))

I want to have the below output to be in my final dataframe which i want to use for renaming the original file.
 Run        source_name
   <chr>      <chr>
 1 SRR2920466 HSC1
 2 SRR2920477 HSC2
 3 SRR2920478 HSC3
 4 SRR2920506 HSC4
 5 SRR2920507 HSC5
 6 SRR2920532 HSC6
 7 SRR2920531 HSC7
 8 SRR2920469 CMP1
 9 SRR2920470 CMP2
10 SRR2920481 CMP3

Found a solution but I have to apply it to each string in the last column
i1 <- !sapply(df1, is.numeric)
df1[i1] <- lapply(df1[i1], function(x) factor(replace(as.character(x),
                                                      x == "HSC", paste0(x[x=="HSC"], seq_len(sum(x == "HSC"))))))

Right now I am able to match repeated string in the column and add increment.How to apply it for all the other string .

We need to group by source_name and add row id, here is using data.table package:
library(data.table)

# convert to data.table
setDT(df1)

# group by source_name, and add rowid
df1[, myNewColumn := paste0(source_name, rowid(source_name))]

df1
#            Run source_name myNewColumn
#  1: SRR2920466         HSC        HSC1
#  2: SRR2920477         HSC        HSC2
#  3: SRR2920478         HSC        HSC3
#  4: SRR2920506         HSC        HSC4
#  5: SRR2920507         HSC        HSC5
#  6: SRR2920532         HSC        HSC6
#  7: SRR2920531         HSC        HSC7
#  8: SRR2920469         CMP        CMP1
#  9: SRR2920470         CMP        CMP2
# 10: SRR2920481         CMP        CMP3
# ...

See here for more solutions.

EDIT: Wheres the OP accepted this answer, @zx8754's answer below, as well as the ones in the link they gave, are much more elegant and straightforward.

Here is a work-around with table(), lapply() and seq_len() (I have basically used your lapply() call on the output of a table() call):
my_counts <- table(df1$source_name)

> my_counts

Blast   CMP   GMP   HSC   LSC  Mono  pHSC
   23     8     7     7     8     6    19

df1 <- df1[order(df1$source_name),]

> head(df1)
          Run source_name
37 SRR2920546       Blast
38 SRR2920547       Blast
39 SRR2920548       Blast
40 SRR2920549       Blast
41 SRR2920550       Blast
42 SRR2920551       Blast

> names(unlist(lapply(my_counts, seq_len)))
 [1] "Blast1"  "Blast2"  "Blast3"  "Blast4"  "Blast5"  "Blast6"  "Blast7"  "Blast8"  "Blast9"  "Blast10" "Blast11" "Blast12" "Blast13"
[14] "Blast14" "Blast15" "Blast16" "Blast17" "Blast18" "Blast19" "Blast20" "Blast21" "Blast22" "Blast23" "CMP1"    "CMP2"    "CMP3"
[27] "CMP4"    "CMP5"    "CMP6"    "CMP7"    "CMP8"    "GMP1"    "GMP2"    "GMP3"    "GMP4"    "GMP5"    "GMP6"    "GMP7"    "HSC1"
[40] "HSC2"    "HSC3"    "HSC4"    "HSC5"    "HSC6"    "HSC7"    "LSC1"    "LSC2"    "LSC3"    "LSC4"    "LSC5"    "LSC6"    "LSC7"
[53] "LSC8"    "Mono1"   "Mono2"   "Mono3"   "Mono4"   "Mono5"   "Mono6"   "pHSC1"   "pHSC2"   "pHSC3"   "pHSC4"   "pHSC5"   "pHSC6"
[66] "pHSC7"   "pHSC8"   "pHSC9"   "pHSC10"  "pHSC11"  "pHSC12"  "pHSC13"  "pHSC14"  "pHSC15"  "pHSC16"  "pHSC17"  "pHSC18"  "pHSC19"

df1$new_col <- names(unlist(lapply(my_counts, seq_len)))

> head(df1)
          Run source_name new_col
37 SRR2920546       Blast  Blast1
38 SRR2920547       Blast  Blast2
39 SRR2920548       Blast  Blast3
40 SRR2920549       Blast  Blast4
41 SRR2920550       Blast  Blast5
42 SRR2920551       Blast  Blast6

You can store the initial row order somewhere to have resulting data frame in the same order the initial one.

While my answer solves your problem, the data.table solution suggested by @zx8754 is much more elegant and should be the accepted answer in my opinion. I have up-voted that one anyway.