text:
Python module bamtags
Does anyone know this python module, bamtags? I'm using a pipeline (Dovetail Genomics, HiRise) that imports it, but it seems it doesn't exist.
It's  impossible to install it and there is no information on the internet. Does anyone know it?
Additional information:
Hirise is a 3 years old pipeline that uses Chromossome Capture (HiC) data to identify misassembled regions in a genome assembly and then rescaffold them. By integrating 2 different sequencing methods, the confidence of the genome assembly is thus increased.
The workflow is free (we only paid for the data generation). Unfortunately I have no idea what the module does, as I find no information on internet and the scripts are poorly documented. f python-reading is correct, it looks like that is an internally defined module (class) which is in here. So basically, it's there. I am not a python expert so maybe others more experience with python development will confirm or deny this.
The only function that is used is the following: BamTags.mate_mapq(aln), which is inside a loop in aln: for aln in sam.fetch(reference=c1)
Thanks,
Ricardo

Answer from @fridaymeetssunday, converted from comments:
If python-reading is correct, it looks like that is an internally defined module (class) which is in here. So basically, it's there.
Response from poster:
Thank you, guys, so helpful! Yes, that's it. I wasn't aware of the simplicity of import function. I thought it only worked with builtin or pip/conda installed modules. I just had to add the scripts folder to my PATH variable or have the functions in the working directory.


text:
Colinearity analysis with respect to gene orientation
I was reading the manual of MCScanX and found that besides a all-vs-all protein blast file the tool takes a gtf file with four columns only: scaffold name, protein name, from, to.
As the tool is not taking the strand column, it appears to me that probably the most popular tool for genome collinearity analysis is disregarding orientation of genes. Do I get it right? Is there an analogical tool that would perform the collinearity analysis with respect to orientation of genes?
-- edit --
Now I am almost entirely sure that it MCScanX is not taking orientation into account since I found collinear blocks where some of the gene homologs were on the same strand and some of them two different strands (I do hesitate to call them co-"linear" now). So the updated question would be: Is there a good MCScanX alternative that would take orientation of genes into account?
I am not sure if I will improve the question with an example, but i will try anyway. Let's say that A B C D was detected as collinear with A' B' C' D'. The orientation of genes is + + - + for A B C D and + - - - for A' B' C' D'. For sure, these are synteny blocks. The order of genes is the same, but are they collinear? How can I search just for blocks that would be same (or completely revers) orientation?

Sounds like you could just do logic on the from, to columns if they correctly have the from > to if the gene is in the reverse orientation as the reference.

Theoretically it would be possible to encode the orientation information like this, but that would be against gtf specification that has the leftmost and the rightmost positions.

Couldn't that be taken care after the "co-linear" analysis? In your example what would you consider co-linear (B and D with their pairs B' and D')? Sorry it is not my fields but seems interesting

@Llopis Yep, I have ended up in post processing the MCScanX output - I need to post an answer. In the example I made there are two different collinear blocks A+ C- C'- A'+ (reverse duplication) and B+ D+ D-' B'- (a palindrome). Assuming that we allow gaps in between genes.

source: Stack Exchange

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

sort: votes


text:
3D PCA group labelling
I would like to make a 3D PCA but not sure how to label group wise which i can do for 2D PCA
library(ggplot2)
test<- read.csv("NEW_RBP/TF_DISEASE.txt",header = T,row.names = 1,sep = '\t')

Sample = c(rep("HSC",4),rep("Blast",11),rep("LSC",8))
Stage=c(rep("HSC",4),rep("Blast",11),rep("LSC",8))

df_final=data.frame(Sample,Stage)
df_final

pca_data=prcomp(t(test), center=TRUE, scale=TRUE)
pca_data_perc=round(100*pca_data$sdev^2/sum(pca_data$sdev^2),1)

df_pca_data=data.frame(PC1 = pca_data$x[,1],
                       PC2 = pca_data$x[,2], sample = colnames(test))
myColors <- c("red", "black", "blue2")
sample <- as.vector(df_pca_data$sample)
ggplot(df_pca_data, aes(PC1,PC2, colour = factor(Stage)))+#,label=rownames(t(test)))+
  geom_point(size=30)+
  geom_text(aes(label=sample),color="white",size=6,angle =0,parse = TRUE,  family="Bookman", fontface="bold")+
  labs(x=paste0("PC1 (",pca_data_perc[1],")"), y=paste0("PC2 (",pca_data_perc[2],")"))+
  theme_minimal(base_size=20) +
  theme(axis.text.x=element_text(size=rel(3), angle=90))+
  theme(axis.text.y=element_text(size=rel(3), angle=90))+
  theme(axis.title.x = element_text(colour="grey20",size=55,angle=0,hjust=.5,vjust=0,face="plain"))+
  theme(axis.title.y = element_text(colour="grey20",size=55,angle=90,hjust=.5,vjust=.5,face="plain"))+
  scale_color_manual(values=myColors)

The above code give me the figure as i have attached
For 3D PCA I tried this
scores = as.data.frame(pca_data$x)
head(scores)

plot3d(scores[,1:3], col=c(1:4), size=20, type='p',
       xlim = c(-50,50), ylim=c(-50,50), zlim=c(-50,50))
text3d(scores[,1]+2, scores[,2]+1, scores[,3]+1,
       texts=c(rownames(scores)), cex= 0.7, pos=3)

I get something like this as in case of 2D PCA i can label lets say all of the HSC sample as black and so on..Im not sure how to do that for 3D as I tried using factor i couldn;t do
Any help or suggestion would be highly appreciated

You can easily color 3D pca plots in R based on the code given below:
library("scatterplot3d")
colors <- c("#999999", "#E69F00", "#56B4E9") # Number of color according to the number of groups
colors <- colors[as.numeric(iris$Species)] # you can put here the column containing the name of population or sample etc.
pca1 <- prcomp(iris[, -5])
s3d <-scatterplot3d(pca1$x[, 1], pca1$x[, 2],pca1$x[, 3],xlab="Sepal.length",ylab="Sepal.width", zlab="Petal.length", pch = 16,color=colors)
legend("right", legend = levels(iris$Species),
  col =  c("#999999", "#E69F00", "#56B4E9"), pch = 16)

The generated graph is given below:

For further details about how to manipulate plots in 3D have a look at this site.
Edited Response
Using the data given below I plotted a 3D pca that may help solve your problem,
"Col1" "Col2" "Col3" "Col4" "colend"
"H1" 5.1 3.5 1.4 0.2 "HSC"
"H2" 4.9 3 1.4 0.2 "HSC"
"H3" 4.7 3.2 1.3 0.2 "HSC"
"H4" 4.6 3.1 1.5 0.2 "HSC"
"H5" 5 3.6 1.4 0.2 "HSC"
"H6" 5.4 3.9 1.7 0.4 "HSC"
"H7" 4.6 3.4 1.4 0.3 "HSC"
"H8" 5 3.4 1.5 0.2 "HSC"
"H9" 4.4 2.9 1.4 0.2 "HSC"
"H10" 4.9 3.1 1.5 0.1 "HSC"
"H11" 5.4 3.7 1.5 0.2 "HSC"
"H12" 4.8 3.4 1.6 0.2 "HSC"
"H13" 4.8 3 1.4 0.1 "HSC"
"H14" 4.3 3 1.1 0.1 "HSC"
"H15" 5.8 4 1.2 0.2 "HSC"
"B1" 5.7 4.4 1.5 0.4 "blast"
"B2" 5.4 3.9 1.3 0.4 "blast"
"B3" 5.1 3.5 1.4 0.3 "blast"
"B4" 5.7 3.8 1.7 0.3 "blast"
"B5" 5.1 3.8 1.5 0.3 "blast"
"B6" 5.4 3.4 1.7 0.2 "blast"
"B7" 5.1 3.7 1.5 0.4 "blast"
"B8" 4.6 3.6 1 0.2 "blast"
"B9" 5.1 3.3 1.7 0.5 "blast"
"B10" 4.8 3.4 1.9 0.2 "blast"
"B11" 5 3 1.6 0.2 "blast"
"B12" 5 3.4 1.6 0.4 "blast"
"B13" 5.2 3.5 1.5 0.2 "blast"
"B14" 5.2 3.4 1.4 0.2 "blast"
"B15" 4.7 3.2 1.6 0.2 "blast"
"B16" 4.8 3.1 1.6 0.2 "blast"
"B17" 5.4 3.4 1.5 0.4 "blast"
"B18" 5.2 4.1 1.5 0.1 "blast"
"B19" 5.5 4.2 1.4 0.2 "blast"
"B20" 4.9 3.1 1.5 0.2 "blast"
"L1" 5 3.2 1.2 0.2 "LSC"
"L2" 5.5 3.5 1.3 0.2 "LSC"
"L3" 4.9 3.6 1.4 0.1 "LSC"
"L4" 4.4 3 1.3 0.2 "LSC"
"L5" 5.1 3.4 1.5 0.2 "LSC"
"L6" 5 3.5 1.3 0.3 "LSC"
"L7" 4.5 2.3 1.3 0.3 "LSC"
"L8" 4.4 3.2 1.3 0.2 "LSC"
"L9" 5 3.5 1.6 0.6 "LSC"
"L10" 5.1 3.8 1.9 0.4 "LSC"
"L11" 4.8 3 1.4 0.3 "LSC"
"L12" 5.1 3.8 1.6 0.2 "LSC"
"L13" 4.6 3.2 1.4 0.2 "LSC"
"L14" 5.3 3.7 1.5 0.2 "LSC"
"L15" 5 3.3 1.4 0.2 "LSC"

The code is as follows:
library("scatterplot3d")
colors <- c("#999999", "#E69F00", "#56B4E9") # Number of color according to the number of groups
colors <- colors[as.numeric(data.c1$colend)] # you can put here the column containing the name of population or sample etc.
pca1 <- prcomp(data.c1[, -5]) # PCA on columns except the last column
s3d<-scatterplot3d(pca1$x[, 1], pca1$x[, 2],pca1$x[, 3],grid=TRUE,xlab="PC1",ylab="PC2", zlab="PC3", pch = 16,color=colors)
legend("right",legend = levels(data.c1$colend),col =  c("#999999", "#E69F00", "#56B4E9"), pch = 16,inset=-0.04,bty="n")
text(s3d$xyz.convert(pca1$x[, 1:3]), labels = rownames(data.c1),cex= 0.7, col = "black",pos=2.5)

In this way we get following plot:
Hope this helps!!!

welll i have seen that ,but my data is as such if i can describe ,im trying to plot the prcomp object x which contains my samples in the rows and PC in the columns...now sure how to group it as i did it for 2D pca..

So you mean to say that your problem has been solved?

no no ..i was saying i can;t do what i saw in the example ...if you can tell me how to proceed ?i would be glad

Is it possible for you  to share your data?

of course...i will update my question with the data link of the data https://file.io/q8UeZ9

To make it more clear, my H labeled samples are HSC ,B as Blast and L are LSC samples so I want to group them as HSC ,blast and LSC respectively

The link is not working.

okay sorry about it i thought you would see it fast so it was a temp storage..here is the link https://drive.google.com/file/d/19niHtbWL3oCPh8a9IdzzoXIzyWgdf7ao/view?usp=sharing

@krushnachChandra I have edited my response, by arranging your data in the same way you can achieve your desired results.

wow looks cool so the idea is to arange the data and make the the last column for labeling, did you transpose the data as i see my sample name are in rows now..what are the columns?

yes I transpose the data, but it is not your data it is iris dataset from R base and column name does not matter here except for the last column i.e colend but what matters is row name.

yes now i get it but im curious when i transpose this is my dimension   23 1328 ,23 rows and 1328 columns..can you suggest me how can i add a new column i.e colend one way i can open the transposed file and add  a column in R like make group for each cell type replicate ,can you suggest other way? lets say if i have 10000 columns after transposing i dont think that can be opened excel..!!

check documentation of cbind() and rep() functions in R

my goodness i miss the simple way ...thank you...i was thinking in a complicated way....thanks a lot..

finally i made did it..as you suggested...

@AmmarSabirCheema i would like to have a look at this question ...https://stackoverflow.com/questions/72958698/adding-group-information-to-3d-plot-in-factoextra

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


text:
Changing a wide range of colours to a limited gradient
I returned a FeaturePlot from Seurat to ggplot. My plot has a weird range of colours as below

I produced this plot by this code
> head(mat[1:4,1:4])
             s1.1        s1.2 s1.3       s1.4
DDB_G0267178    0 0.009263254    0 0.01286397
DDB_G0267180    0 0.000000000    0 0.00000000
DDB_G0267182    0 0.000000000    0 0.03810585
DDB_G0267184    0 0.000000000    0 0.00000000
>

I have converted expression matrix to a binary matrix by 2 as a threshold

mat[mat < 2] <- 0
mat[mat > 2] <- 1

> head(exp[1:4,1:4])
             s1.1 s1.2 s1.3 s1.4
DDB_G0267382    0    0    0    1
DDB_G0267438    0    0    0    1
DDB_G0267466    0    0    0    0
DDB_G0267476    0    0    1    0
>
> exp=colSums(exp)
> exp=as.matrix(exp)
> colnames(exp)="value"
> exp=as.data.frame(exp)
> cc <- AddMetaData(object = seurat, metadata = exp)
> cc <- SetAllIdent(object = cc, id = "value")
> TSNEPlot(object = cc, do.return= TRUE)

How I can convert this range to a gradient of colours for example in 8-18, 18-28, 28-38, 38-48 range in blue to yellow please? Something like below

Thank you for any help
Then by ggplot now I scaled my colours but I don't like my clusters as so and I don't know how to retain my clusters as a featureplot by this new color gradient
> head(cc@meta.data)
     nGene    nUMI    orig.ident res.0.7 CELL STAGE GENO dataset stage.nice celltype value
s1.1  4331  373762 SeuratProject       0 s1.1   H16   WT       1        H16        0    34
s1.2  5603 1074639 SeuratProject       0 s1.2   H16   WT       1        H16        0    26
s1.3  2064   49544 SeuratProject       0 s1.3   H16   WT       1        H16        0    27
s1.4  4680  772399 SeuratProject       1 s1.4   H16   WT       1        H16        1    29
s1.5  3876  272356 SeuratProject       1 s1.5   H16   WT       1        H16        1    21
s1.6  2557  122314 SeuratProject       0 s1.6   H16   WT       1        H16        0    31

> ggplot(as.data.frame(cc@meta.data), aes(x = cc@meta.data$CELL, y = cc@meta.data$res.0.7, colour =cc@meta.data$value)) +
+     geom_point(size = 5) +
+     scale_colour_gradient(low = "yellow", high = "blue")

By below code I obtained a tsne in link
> cols <-  scales::seq_gradient_pal(low="beige", high="red", space="Lab")(seq(from=0, to=1,length.out=48))
>
> TSNEPlot(cc, colors.use=cols)

enter link description here
Now I want to know how I could convert this range to a 8-18, 18-28, 28-38, 38-48 colour range as a gradient of blue to yellow?

Can you share the code used to generate the first plot? Essentially what you put in to the call to FeaturePlot().

Yes, Sure; I will add to the main question

Looks like your value column is not treated as an integer but as a factor.

Coloring aside, what are trying to accomplish with >colnames(exp)="value"? Unless I'm missing something you are trying to assign one string to a vector of values. This shouldn't work. You should be getting an error like length of 'dimnames' [2] not equal to array extent. Do you ultimately want a plot of the values in one column of the matrix named exp?

I missed some part, I added that. I have used colSums of binary matrix to be able to add that to seurat metadata; exp=colSums(exp)    Please don't negatively point me, my problem is really difficult to solve

The exp and the mat matrices do not match. Please show a reproducible, minimal example.  Also by converting the exp to a data.frame it might convert your numeric values to factors, which is why you see them with each factor into a different color.

Actually I done that by ggplot, I will add to my question

@Llopis, the data tables ought not match. exp is a binarized version of mat as Feresh Teh states above.
Feresh Teh: I would think you would want rowSums(), that way all of the columns for a given row will be reduced to the sum and you'll have a value (expression?) for each cell. When you add metadata to a Seurat object, it will be in this format--a value for each cell.

@kohlkopf Unless I misunderstood exp is exp <- mat;exp[exp>2] <- 1;exp[exp<2] <- 0 But they don't match the row names between the two heads shown.

Ah, you are pointing out the row names specifically. We are definitely missing some of the steps here, a minimally reproducible example is a good idea.

I previously thought the the column names were genes. Please disregard, colSums() makes sense. The cell names are the s1.1, s1.2 etc. Working on an answer.

If you would like to color discrete intervals on a gradient as opposed to having a continuous gradient (like your second plot), use this approach.
It is similar to the approach in the answer I posted with the continuous scale, but we simply break up the continuous scale in to intervals and color them by these intervals.
#generate values for testing purposes, one value for each cell, add to object
value <- sample(seq(from=8, to=46, by=1), size = length(rownames(exp@meta.data)), replace=TRUE)
exp <- AddMetaData(object=exp, metadata=value, col.name="value")

#encode the continuous values as factors, determined by the interval they fall in to
value_breaks <- cut(exp@meta.data$value, breaks = c(8,18,28,38,46), include.lowest=TRUE, right=FALSE)

#name the breaks by cell so they can be added by AddMetaData(), add them
names(value_breaks) <- rownames(exp@meta.data)
exp <- AddMetaData(object=exp, metadata=value_breaks, col.name="value_breaks")

exp <-  SetIdent(exp, ident.use=exp@meta.data$value_breaks)

#create a gradient between specified intervals
#length.out is the number of different colors, number of factor levels
cols <- scales::seq_gradient_pal(low="blue", high="yellow", space="Lab")(seq(from=0, to=1,length.out=4))

TSNEPlot(exp, colors.use=cols)

The label in the legend are classic mathematical notation for intervals. You can add your own labels if you wish in the cut() function.

Thanks a lot, you solved my long lasting question, I produced this final picture https://image.ibb.co/dsZzsf/Rplot22.png  here I am plotting 300 cell types markers hopefully to gene 2 clusters yellow and another one blue, but as you are seeing the threshold I have selected for converting mat (expression values) of these genes to a binary matrix has not been good because my tsne does not look brilliant. If you where me, what threshold you would selected to converting your expression values to binary? I have used the average of expression of these genes in clusters as threshold. for example 2

That deserves a post of its own! That exact question was asked on here a few weeks back, but I think it was removed. There were at least four different "answers." It is a difficult questions with all answers (that I've seen) with their own limitations.

Thank you, in CellRouter R package, developer uses a centring on data before plotting but I am not able to understand what he is doing. But his centring is very good https://github.com/edroaldo/cellrouter/blob/master/CellRouter_Class.R  in plotDRExpression  function

TSNEPlot()
TSNEPlot() will always treat your variables as discrete. My approach is to manually generate a gradient with unique colors for each factor level and pass it to the cols.use argument in TSNEPlot().
#generate values for testing purposes, one value for each cell
value <- sample(seq(from=8, to=48, by=1), size = length(rownames(unfiltered_cca@meta.data)), replace=TRUE)
names(value) <- rownames(exp@meta.data)
exp <- AddMetaData(object=exp, metadata=value, col.name="value")
exp <-  SetIdent(exp, ident.use=exp@meta.data$value)
#create a gradient between specified colors, multiply by sequence to get appropriate length
#length.out is the number of different colors, number of factor levels
cols <-  scales::seq_gradient_pal(low="beige", high="red", space="Lab")(seq(from=0, to=1,length.out=48))

TSNEPlot(exp, colors.use=cols)

I would present this with no legend:
TSNEPlot(exp, colors.use=cols, no.legend = TRUE)

FeaturePlot()
You can also simply use FeaturePlot() instead of TSNEPlot() to visualize the gradient. Using the same data as above:
FeaturePlot(object = exp, features.plot = "value", reduction.use = "tsne", no.legend = FALSE, cols.use = c("beige", "red"))

You ask for a continuous scale, but this is not what is shown in your second plot.  You have a different color for each discrete range of values. I assume here that you want a continuous gradient. To achieve something like the second plot you provided requires a different approach.

Sorry, your solution returned colourless featureplot

TSNEPlot() returns a colourless plot. I have used cc@meta.data$value as my value and in col I used cols <- scales::seq_gradient_pal(low="beige", high="red", space="Lab")(seq(from=10, to=48, length.out=34)), however returned a colourless lot

seq(from=0, to=1, length.out=34) is essential. You need a sequence of decimals for the color gradient. The length.out is correct.

Actually like my second plot I would like to have a colour range for example 8-18, 18-28, 28-38, 38-48 in gradient of colours but I don't know how to do that

Thanks a lot, now I have a colourful tsne plot, sorry do you know how to convert this wide range of colours to 4 ranges 8-18, 18-28, 28-38, 38-48 as a yellow to blue gradient like my second plot?

Sorry how I could convert cols <-  scales::seq_gradient_pal(low="beige", high="red", space="Lab")(seq(from=0, to=1,length.out=48))
to theme_bw() +
#scale_colour_gradientn(name = "Expression", colours=rev(rainbow(4))) +
scale_colour_gradientn("Relative expression", colours=c("midnightblue","dodgerblue3","white", "darkorange2", "yellow"))


text:
How to get the corrected matrix after SVA batch effect correction
I ran SVA to remove batch effects for my bulk RNAseq experiments, but now I need to somehow correct my data matrix in order to run pca, mds. I am using DESeq2 for the analysis. Here is the code that I got now:
dds <- estimateSizeFactors(dds)
dat = counts(dds, normalized = TRUE)
idx = rowMeans(dat) > 1
dat = dat[idx,]
mod = model.matrix(~Group, colData(dds))
mod0 = model.matrix(~1, colData(dds))

nsv = num.sv(dat, mod)
svseq = svaseq(dat, mod, mod0, n.sv = nsv)

dds_sva = dds
dds_sva$SV1 = svseq$sv[,1]
design(dds_sva) = ~ SV1 + Group
dds_sva = DESeq(dds_sva)

So, I am wondering now, how to use svseq object to plot, say, the corrected version of MDS plot.

If you want to plot the "corrected" expression, you will need to remove the variation introduced by these surrogate variables. Removing the expression affected can introduce some bias too and it is usually not recommended (despite comBat doing so). You should apply linear algebra, you can look at here is an example how to do it:
library("corpcor")
s <- corpcor::fast.svd(t(scale(t(dat), center = TRUE, scale = TRUE)))
pcSds <- s$d
pcSds[1] <- 0
svdexp <- s$u %*% diag(pcSds) %*% t(s$v)
colnames(svdexp) <- colnames(dat)

But if your surrogate variables do have a big enough influence on your MDS plot to be observable (IMHO) there would be too much batch effect to work with the data.

But batch effect often appears as grouping of points, and that is how often I can see it. So, I am a bit confused regarding it. I would expect this grouping to disappear and get a different clustering of points. Are you talking only about MDS or about pca too?

Well, without knowing the data is difficult to know if there is really a batch effect or other variables that could be adjusted for. I'm talking about both PCA and MDS, usually there isn't such a big difference between them (unless you use some other metric for distance in the MDS that is not the eucledian distance).


text:
Delete all 4 lines of a fastq read from a fastq file using read ID
I have the following error when running bowtie2:

Error: Read HWI-D00466:116:CC62WANXX:3:1102:7363:63646 1:N:0:GCACACG
  has more read characters than quality values.

I now want to remove all 4 lines of this specific read from the fastq file.
How can I use awk or sed to do this?

Run tail on your FASTQ file. If a job aborted while writing to disk, this read will be the last line. This error message likely means the file is truncated and missing reads. That read is not as much of a concern as the rest of the data. It’s best to run your data processing again to make sure the file is complete.

Don't do it: Your FASTQ file is either malformed or your FASTQ record spans more than four lines, which is allowed in FASTQ.  For a detailed description of what can go wrong in FASTQ parsing see for example http://biopython.org/DIST/docs/api/Bio.SeqIO.QualityIO-module.html#FastqGeneralIterator.
If the FASTQ is malformed, then you should really ask yourself how this happened in the first place and fix the source of the problem. If the record is valid FASTQ, then I suggest parsing the read with for example FastqGeneralIterator and dumping the parsed result back to FASTQ in a 4-line-per-record form.

If you are 100% sure the read only has 4 lines (they can have more), you can use this sed command:
sed -i.bak '/^@HWI-D00466:116:CC62WANXX:3:1102:7363:63646 1:N:0:GCACACG/,+3d'

The -i.bak makes sed modify the original file and create a backup copy with the same name and the extension .bak. The command just means "delete the line matching the pattern and the next three lines".

Maybe the "@" should be added to the pattern. What would happen if the read name was also on the "+" line?

When I can I like to do file processing line-by-line in the spirit of UNIX tools. You can read 4 lines from a Fastq file into 4 tab-separated values on a single line using paste, and then use grep to filter out the record in question. (Then you just have to turn the tabs back into line breaks.)
paste - - - - < reads.fastq \
    | grep -v 'HWI-D00466:116:CC62WANXX:3:1102:7363:63646' \
    | tr '\t' '\n' \
    > reads-fixed.fastq

However, it's likely that other reads in your Fastq file are corrupted as well, in which case it's probably better to write a script in Python or some other language that discards all reads whose read length doesn't match the quality string length.
But of course the most important question is: how did the data become corrupted in the first place? Where did this file come from, and what pre-processing steps were performed on it? Is getting rid of this read (or all reads with mismatched sequence/quality length) going to solve the problem for you? These aren't questions that we can answer for you, but probably deserve good answers. :-)

Here's an awk one-liner that works as intended, unlike the bioawk answer.
This removes all instances of fastq records in which the length of the seq and qual fields do not match. Again, append | gzip > filtered.fastq.gz to gzip the output.
zcat bad_file.fq.gz | awk '{ pos = NR%4; if ( pos == 1) {h=$0} else if (pos == 2) \
  {s=$0} else if (pos==3){c=$0} else if ( pos==0 ) \
  {q=$0; if ( length(q) == length(s) ){printf("%s\n%s\n%s\n%s\n", h, s, c, q)}}}'

Delete the \ characters at the end of the first two lines to make it easy to paste into your terminal.
Here is a shorter version of the same thing thanks to @terdon
zcat file.fastq.gz | awk '{ record[++k]=$0; if(NR%4==2){slen=length($0)} else if(NR%4==0){if(length($0) == slen){ for(i in record){print record[i]}} k=0; }}'

Heh, I was just writing one too :) Here's what I came up with, it's a bit shorter than yours: zcat file.fastq.gz | awk '{ record[++k]=$0; if(NR%4==2){slen=length($0)} else if(NR%4==0){if(length($0) == slen){ for(i in record){print record[i]}} k=0; }}' feel free to add it to your answer if you want to.

To generalize this to remove all reads that do not have the same length of read and sequence quality, you can use this bioawk one-liner:
bioawk -cfastx '{ if (length($seq) == length($qual) ){printf("@%s %s\n%s\n+\n%s\n", $name, $comment, $seq, $qual)} }' my_fastq.gz

This strips anything in the "fastq comment field", +, but that almost never has anything there anyway.
Append | gzip > filtered.fastq.gz to the above command to gzip the output.

This doesn't seem to work for me (Arch Linux, bioawk version 20110810). Even after correcting the syntax error (corrected in your answer), it seems to exit after the first record it finds where the lengths are not the same. Can you confirm it works on your machine?

Thanks for the edit! Yep, I just confirmed that it works for me on all records. I'm using Ubuntu 16.04. >>bioawk --version returns awk version 20110810. It looks like the last time that I compiled bioawk was from the github repo in 2016.

That's really weird. I have the same version. Can you try a case where the record with different lengths is the second one in the file? Does it not just print the first record there?

I got the same error as you - I made the second entry in the fastq file have a different length of seq and qual fields and bioawk just exited before printing that record. I wonder if this is intended behavior? I have to add that in my above comment, I just tried the bioawk one-liner on a normal fastq file without any length mismatches. Oops!

I don't see how it can be intended. Sounds like a bug in bioawk.


text:
Total reads aligning to each reference within a bam file
I have two PCR amplicons that have been multiplexed and sequenced using the nanopore minion.
I have aligned the fastq reads using minimap2 with a reference file containing both amplicon sequences and generated a bam file that I have viewed using IGV.
I am looking for a way to generate some simple summary statistics.
In particular, is there a way to extract the total number of fastq reads aligning to each amplicon reference from the bam file?

The quick way to get the number of alignments on each reference is
samtools idxstats my_bam.bam

Number of reads on each reference is column 3. Although, as has been pointed out, this will give you the total number of alignments per reference, not the total number of reads (each read might give rise to more than one alignment). That said I do tend to us this as generally I'm after a rough approximation, rather than an accurate number.
In theory, only one alignment for each read should be marked as primary, so the following should give you what you need quickly and at low memory usage:
samtools view -bF 2304 my_bam.bam > primary_only.bam
samtools index primary_only.bam
samtools idxstats primary_only.bam

This one-liner below will work better for long reads than samtools flagstat in that it only counts the primary alignment for each read and samtools flagstat doesn't seem to calculate some stats for long reads. I have never seen samtools flagstat output stats on a per-reference basis, but am curious if so!
This answer filters out secondary and supplementary alignments for your reads (-F 2304) that have some alignment to both amplicon reference and just keep the best one. This might give a more accurate idea of how many reads of each amplicon are in the library.
samtools view -F 2304 myfile.bam | awk -F $'\t' '{a[$1, $3]++} END{for (i in a) {split (i, sep, SUBSEP); print sep[1], sep[2], a[i]}}' | uniq | awk '{print($2)}' | uniq -c | sort -k1 -nr

adapted from this.

This can be done using
   samtools flagstat your_bam_file

I think that flagstat counts the number of alignments in the bam file.


text:
Using Beagle 4.1 for phasing and IBD - chromosomes
I am trying to use beagle 4.1 for phasing and IBD analyses. I was wondering if I need to run the analysis by chromosome or I can put them together in the same vcf file.

Well, I am new here, do you mind to give me some feedback on why the -1 on my question? Not trying to argue, just trying to improve.

I'm not sure who downvoted your question, but maybe because you don't have any info about what you have done so far or what you aim to achieve? In any case I think the question and the answer you need are clear enough, so +1. If you find out the answer to your question by some other means we would be happy if you submit an answer to your question.

Thank you for your explanation. I believe it should be clear enough for anyone familiar with the program. I ended up contacting one of the authors, I will post his answer here.

OK, I ended up contacting one of the authors of the program, who kindly answered to my question:

Either approach should work,

[i.e. running the analysis by chromosome vs joining the data together]

but I suggest dividing by chromosome because that approach has had
  more testing.

I hope this might help other users, I saw that several questions regarding Beagle went unanswered on biostars and seqanswers. I guess it is a pretty specialized topic.
Best,
Max

I don't see many questions about phasing and IBD detection on bioinformatcis SE. I think that's just because it's still a small site and there aren't that many people from that field on it yet. I wish there were more questions!

I would also add that splitting your data into chromosomes allows you to run each chromosome in parallel (presuming you have some kind of server you can send each run to), vastly reducing the amount of time your analysis is going to take. I've not run all 22 chromosomes at once, but I can imagine it will take an enormous amount of time if you are using even a small number of markers and references.

This would be better served as a comment to the answer since it does not stand as a answer itself


text:
GO Term heatmap plot in terms of P value or fold enrichment
I'm clustering genes in terms of expression after clustering them. I'm taking out clusters and trying to find out what kind of GO terms are coming up. I came across this figure from this paper, I want to do as the figure B:
Here you can find the file I intend to use to make something like the above plot.
I am using this the code from a question in Biostar:
DAVIDfile <- "DAVID.txt"
DAVID <- read.table(DAVIDfile, sep="\t", header=TRUE)
colnames(DAVID)
class(DAVID)

enrichBcutoff <- 0.05
DAVID <- subset(DAVID, Benjamini< enrichBcutoff)
DAVID <- DAVID[,c(1,2,3,5,6,10)]

head(DAVID)

A <- data.frame(sigGeneList)

annGSEA <- data.frame(row.names=A)
for (j in 1:length(sigGeneList))
{
  gene <- sigGeneList[j]
  pattern <- paste("^", gene, ", |, ", gene, "$| ", gene, ",", sep="")
    for (k in 1:nrow(DAVID))
  {
    if (any(grepl(pattern, DAVID$Genes[k])))
    {
      annGSEA[j,k] <- 1
    }
    else
    {
      annGSEA[j,k] <- 0
    }
  }
}
colnames(annGSEA) <- DAVID[,2]

annGSEA <- annGSEA[,apply(annGSEA, 2, mean)!=0]

annGSEA <- annGSEA[apply(annGSEA, 1, mean)!=0,]

require(ComplexHeatmap)
require(circlize)

topTable <- topTable[which(topTable$Gene %in% rownames(annGSEA)),]
topTable <- topTable[match(rownames(annGSEA), topTable$Gene),]

geneLab=10
termLab=8

#Create heatmap annotations
dfMinusLog10FDRGenes <- data.frame(-log10(topTable[which(topTable[,1] %in% rownames(annGSEA)),"padj"]))
dfMinusLog10FDRGenes[dfMinusLog10FDRGenes=="Inf"] <- 0
dfFoldChangeGenes <- data.frame(topTable[which(topTable[,1] %in% rownames(annGSEA)),"log2FoldChange"])
dfGeneAnno <- data.frame(dfMinusLog10FDRGenes, dfFoldChangeGenes)
colnames(dfGeneAnno) <- c("DEG\nsignificance\nscore", "Regulation")
dfGeneAnno[,2] <- ifelse(dfGeneAnno[,2]>0, "Up-regulated", "Down-regulated")
colours <- list("Regulation"=c("Up-regulated"="royalblue", "Down-regulated"="yellow"))
haGenes <- rowAnnotation(df=dfGeneAnno, col=colours, width=unit(1,"cm"))

dfMinusLog10BenjaminiTerms <- data.frame(-log10(read.table(DAVIDfile, sep="\t", header=TRUE)[which(read.table(DAVIDfile, sep="\t", header=TRUE)$Term %in% colnames(annGSEA)),"Benjamini"]))
colnames(dfMinusLog10BenjaminiTerms) <- "GO Term\nsignificance\nscore"
haTerms <- HeatmapAnnotation(df=dfMinusLog10BenjaminiTerms,
                             colname=anno_text(colnames(annGSEA), rot=40, just="right", offset=unit(1,"npc")-unit(2,"mm"), gp=gpar(fontsize=termLab)),
                             annotation_height=unit.c(unit(1, "cm"), unit(8, "cm")))

pdf("GO.pdf", width=7, height=12)
hmapGSEA <- Heatmap(annGSEA,

                    name="My enrichment",

                    split=dfGeneAnno[,2],

                    col=c("0"="white", "1"="forestgreen"),

                    rect_gp=gpar(col="grey85"),

                    cluster_rows=T,
                    show_row_dend=T,
                    row_title="Statistically-significant genes",
                    row_title_side="left",
                    row_title_gp=gpar(fontsize=12, fontface="bold"),
                    row_title_rot=0,
                    show_row_names=TRUE,
                    row_names_gp=gpar(fontsize=geneLab, fontface="bold"),
                    row_names_side="left",
                    row_names_max_width=unit(15, "cm"),
                    row_dend_width=unit(10,"mm"),

                    cluster_columns=T,
                    show_column_dend=T,
                    column_title="Enriched terms",
                    column_title_side="top",
                    column_title_gp=gpar(fontsize=12, fontface="bold"),
                    column_title_rot=0,
                    show_column_names=FALSE,
                    #column_names_gp=gpar(fontsize=termLab, fontface="bold"),
                    #column_names_max_height=unit(15, "cm"),

                    show_heatmap_legend=FALSE,

                    #width=unit(12.5, "cm"),

                    clustering_distance_columns="euclidean",
                    clustering_method_columns="ward.D2",
                    clustering_distance_rows="euclidean",
                    clustering_method_rows="ward.D2",

                    bottom_annotation=haTerms)

draw(hmapGSEA + haGenes, heatmap_legend_side="right", annotation_legend_side="right")
dev.off()

However there is an error message after annGSEA:
dput(annGSEA)
structure(list(), .Names = character(0), row.names = "c(\"LRP1\", \"ELAVL4\", \"ENOX1\", \"TRIP6\", \"KHDRBS3\", \"TDRD9\", \"THOC6\", \"CPSF1\", \"HSPB1\", \"SF3A2\", \"STXBP1\", \"DUS3L\", \"INTS1\", \"DZIP1L\", \"HIST1H1B\", \"TRMT2A\", \"RPUSD1\", \"MRPS2\", \"EXOSC5\", \"TDRD1\", \"PURG\", \"NOL12\", \"TRMT61A\", \"XAB2\", \"TRMT1\", \"PTRF\", \"POLR2I\", \"MAZ\", \"TBL3\", \"BOP1\", \"RBM10\", \"TSR3\", \"SART1\", \"CCDC86\", \"SAMD4A\", \"RBM11\", \"RPS5\", \"HNRNPLL\", \"POLRMT\", \"NSUN7\", \"PABPC5\", \"DDX54\", \"PRRC2B\", \"RPS15\", \"MYBBP1A\", \"DUS1L\", \"ARHGEF28\", \"RPL28\", \"ZNF598\", \"ZMAT5\", \"RPL36\", \"RPL8\", \"DHX37\", \"ZC3H7B\", \n\"RRP9\", \"PQBP1\", \"ZNF385A\", \"PUS1\", \"DDX49\", \"TCOF1\", \"BICC1\", \"GLTSCR2\", \"SCAF1\", \"MRPL4\", \"FSCN1\", \"EDF1\", \"NOVA1\", \"PCBP3\", \"FBL\", \"MRPL55\", \"PUF60\", \"NYNRIN\", \"ZFP36L1\", \"RPS9\", \"MKRN3\", \"GTF2F1\", \"EEFSEC\", \"FARSA\", \"TARBP2\", \"RPL35\", \"RPS19\", \"ZC3H3\", \"H1F0\", \"BST2\", \"FASN\", \"ZYX\", \"NOC4L\", \"LONP1\", \"RUVBL2\", \"MRPL38\", \"TRAP1\", \"TUFM\", \"RBM42\", \"TNS1\", \"EDC4\", \"WDR6\", \"PCBP4\", \"PUSL1\", \"EEF2\", \"RPS28\", \"RAVER1\", \"APEH\", \"TSEN54\", \"VARS\", \"C7orf50\", \"GRN\", \"ARHGEF1\", \"SF3B5\", \"RPL18\", \"EIF3G\", \n\"ANG\", \"EIF5AL1\", \"HIST1H1C\", \"RPL13\", \"U2AF1L4\", \"C9orf114\", \"FASTK\", \"TSEN34\", \"SRRT\", \"RPS2\", \"MRPL54\", \"CSRP1\", \"HELZ2\", \"RPLP2\", \"RPL18A\", \"PRRC2A\", \"GTPBP1\", \"NOC2L\", \"RBM19\", \"RNASEH2A\", \"ADARB2\", \"SERPINH1\", \"CCDC9\", \"DXO\", \"GNB2L1\", \"FLNA\", \"RNASEH2C\", \"IPO13\", \"HABP4\", \"C11orf68\", \"QTRT1\", \"DGCR14\", \"RPL3\", \"SETD1A\", \"GEMIN4\", \"NPM2\", \"MRPS26\", \"RRP7A\", \"CCDC124\", \"CACTIN\", \"ISG20\", \"INTS5\", \"SNTB2\", \"RPUSD3\", \"RPL29\", \"PPRC1\", \"PPAN\", \"PRPF40B\", \"RBFOX2\", \"C15orf52\", \"ALYREF\", \"EXOSC4\", \n\"GNL1\", \"SIDT2\", \"REPIN1\", \"CHERP\", \"NOB1\", \"RPLP1\", \"EIF5A2\", \"RRP1\", \"MRPL41\", \"RPL19\", \"SUGP1\", \"NAA38\", \"CDC42EP4\", \"EEF1D\", \"DDX41\", \"PRPF19\", \"CLASRP\", \"PRPF6\", \"ZC3H12B\", \"SNRPB\", \"MEPCE\", \"COL14A1\", \"WBSCR16\", \"LSM7\", \"PELP1\", \"VARS2\", \"TRPT1\", \"RPS10\", \"PTBP1\", \"RPS8\", \"CPSF3L\", \"NOSIP\", \"TCF20\", \"PTCD1\", \"APEX1\", \"C1orf35\", \"MRPL23\", \"TWF2\", \"CTIF\", \"DHX38\", \"MRPS34\", \"SYMPK\", \"U2AF2\", \"QARS\", \"CCAR2\", \"RPS3\", \"RPS16\", \"P4HB\", \"MVP\", \"POLR2J\", \"SKIV2L\", \"GTPBP3\", \"ZNF768\", \"UBA52\", \"SUPT5H\", \n\"POLR2L\", \"POLR2A\", \"EIF3K\", \"R3HCC1\", \"PIH1D1\", \"ATXN2L\", \"OAS3\", \"IPO4\", \"CCDC137\", \"SETD1B\", \"ILF3\", \"TBRG4\", \"PKM\", \"AARS2\", \"FAU\", \"ACTN4\", \"RPS19BP1\", \"ANGEL1\", \"PLD6\", \"MOV10\", \"EIF6\", \"DRG2\", \"AARS\", \"PRR3\", \"AC004381.6\", \"RBPMS\", \"CPSF4\", \"NOL6\", \"RPL13A\", \"NPM3\", \"UBA1\", \"REXO4\", \"RRP12\", \"UPF1\", \"RPL11\", \"DDX43\", \"SNRPA\", \"ISY1\", \"MRPS18A\", \"RBMS1\", \"CSTB\", \"FKBP3\", \"MTIF3\", \"MRPL39\", \"CLK4\", \"SRP54\", \"ERI1\", \"CLK1\", \"IMP3\", \"GRSF1\", \"GTF3A\", \"FASTKD1\", \"DYNC1LI1\", \"PTGES3\", \"TRIT1\", \"EXOSC3\", \n\"IGF2BP3\", \"TRIM71\", \"MCTS1\", \"GEMIN2\", \"SUMO1\", \"NQO1\", \"HSPA8\", \"POLR2G\", \"TAF9\", \"MRPL33\", \"FAM46A\", \"MRPS36\", \"FNDC3B\", \"DDX52\", \"DYNLL1\", \"RBM43\", \"MANF\", \"TDRD7\", \"AHNAK\", \"IFIT2\", \"ZMAT3\", \"APOBEC3H\", \"RTN4\", \"S100A4\", \"PIH1D2\", \"ANXA2\", \"PIWIL3\", \"TMSB4X\", \"PIWIL4\", \"RNASE6\", \"LGALS1\", \"RNASE3\", \"TLR7\", \"MAEL\", \"SAMHD1\", \"RNASE2\", \"C2orf15\", \"SIDT1\", \"HERC5\", \"PSIP1\", \"RPL17\", \"ZC3H12C\", \"SLC16A3\", \"TEP1\", \"HEXIM1\", \"CPEB4\", \"ZC3H12D\", \"CORO1A\", \"ZFP36\", \"KIF1C\", \"WARS\", \"CCRN4L\", \"AFF2\", \n\"TST\", \"RBM47\", \"LGALS3\", \"NXF3\", \"PPARGC1A\", \"TLR8\")", class = "data.frame")
>

I think you forgot the link to the paper you are referencing. You tagged the question with [tag:R], which problem have with which code? Perhaps you want to make the figure B and you don't know how?

yes i m trying to make figure and truly i don;t know how to begin with

You added a lot of code, you already seem how to start! What is the output of this code? (could you upload the hmapGSEA plot?) What didn't you like or is not working as expected?

I took this code from biostar, im stuck at the annGSEA after that when i loop through it doesn't match the pattern of genes would give the error message in the question have a look

Please provide the output for the variable pattern, the error points to there being something wrong with that. Also, it looks like you didn't copy paste the error message correctly. Give us the full message and it'll be easier the pin the issue down.

i should use dput to make it clear perhaps.


text:
Gene not found in Affymetrix expression profiles
I am studying the ABA network in A. thaliana, consisting of HB7, ABI1 and AREB2. The AGI code I was given are, respectively: AT2G46680, AT4G26080 and AT1G45249.
I downloaded the following file in order to convert the array element name to the AGI one, but I am not able to find AT1G45249 in the above file. How can I check whether AT1G45249 has another name? Or even better, how can I check what is the array element name in Affymetrix (GPL198) corresponding to AT1G45249. Thank you.

You can use  org.At.tair.db: Genome wide annotation for Arabidopsis, which is an R package available on Bioconductor to extract gene symbols.


text:
What is the minimum number of individuals per group required for ADMIXTURE analysis?
What is the minimum number of individuals required per group (population/species) for an ADMIXTURE analysis to be appropriate?
Is there a way to determine this number mathematically based on the amount of genetic data I have or other criteria?
Are there any citations I can use to figure this information out?

The following paper provides an in-depth exploration of factors affecting power of admixture analyses under a number of different models. Although the authors state they "provide the open-source R software AdmixPower, to perform power and sample size analysis", in fact it seems the software is not available.
AdmixPower: Statistical Power and Sample Size Estimation for Mapping Genetic Loci in Admixed Populations
Tool availability:
https://research.cchmc.org/mershalab/Tools.html
Site asks for a password, apparently there is no option to sign up.


text:
How are the values of prop.part() and prop.clades() calculated?
Consider the following dataset:
fictional.df <- data.frame(L1 = c(0,0,0,0,0,0,0,0),
                       L2 = c(0,1,0,0,0,1,1,0),
                       L3 = c(1,1,0,1,1,1,1,1),
                       L4=c(0,0,1,1,0,0,0,0))

I

converted this to a phyDat object and then
created a pairwise distance matrix as follows:

    fictional.phydat <- as.phyDat(fictional.df,
                              type="USER",levels=c("1","0"),
                              names=names(fictional.df))
    fictional.hamming <- dist.hamming(fictional.phydat)

From this distance matrix, I then estimated a UPGMA tree:

    fictional.upgma <- upgma(fictional.hamming)

I then created bootstrap datasets:

    set.seed(187)
    fictional.upgma.bs <- bootstrap.phyDat(fictional.phydat, FUN =
    function(xx) upgma(dist.hamming(xx)), bs=100)

I then calculated the proportion of partitions in the bootstrap set:

    upgma.bs.part <- prop.part(fictional.upgma.bs)

So far so good. Here is where I would appreciate some help. When I call the function prop.clades, I do not understand the result:

    prop.clades(fictional.upgma,fictional.upgma.bs)
    [1] 100  NA  71

Question Why does this function return NA when there is evidence for that clade in the set of bootstrap trees?
A second question:
prop.clades(fictional.upgma,part=upgma.bs.part)
[1] 100  49 112

If there are only 100 bootstrap samples, why is the value for the final clade 112?


text:
Error given while trying to index a BAM file with Samtools Index - NO COOR?
I am currently working on my own Metagenomic pipeline, utilizing Bowtie 2 to map. Bowtie 2 outputs a SAM file, which I convert to a .BAM and sort it using Samtools. When I try to utilize Samtools to index my .BAM file it gives me this error:
[E::hts_idx_push] NO_COOR reads not in a single block at the end 1652 -1.

Here is my code:
source /opt/asn/etc/asn-bash-profiles-special/modules.sh
module load samtools/1.3.1

cd ~/gz_files/soapdenovo
r=20
###SAM to BAM and sort

samtools faidx soapdenovo.fa
samtools view -b -T soapdenovo.fa -o sample${r}bowtiemapping.bam sample${r}bowtiemapping.sam
samtools sort -n -T ~/gz_files/soapdenovo/temp -o sample${r}bowtiemapping_sorte$
samtools index sample${r}bowtiemapping_sorted.bam

################LETS MAKE THE ABUNDANCE TABLES###########################$
samtools view -c -f 4  sample${r}bowtiemapping_sorted.bam > sample${r}alignment$

samtools idxstats sample${r}bowtiemapping_sorted.bam >  sample${r}idxstats.txt

#git clone https://github.com/metajinomics/mapping_tools.git
#python mapping_tools/get_count_table.py *.idxstats.txt > contig_counts.tsv
##less contig_counts.tsv

Any help would be appreciated!

Notes on pipeline

samtools import was replaced with samtools view
faidxcorrected
Instructions used are here on github
Suggestion for future trouble shooting is given here  by @winni2k

In you code, which is the line that is causing trouble? I cannot see any samtools index for indexing the bam, as you say.

BTW, don't search for "null". This is a name for something that does not exists. Check if all your records have a non-empty sequence.

More comments... samtools import has been deprecated for 8 years. Use samtools view instead.

In your usage of samtools view, you didn't add neither the parameter to tell samtools the input is a sam file nor the parameter to tell samtools the output should be a bam file. By default, you get the other direction.

Finally, the second samtools faidx makes no sense: there are no fasta files to index!

I was following a tutorial from the 2018 Edamame course (https://github.com/edamame-course/2018-Tutorials/wiki/Schedule-EDAMAME-2018)

I am new to bioinformatics so any tips would be helpful. This is the line that is giving me issues:
samtools faidx sample${r}bowtiemapping_sorted.bam sample${r}bowtiemapping_sorted.bam.bai

I made the recommended changes above. I am now getting this error: [E::hts_idx_push] NO_COOR reads not in a single block at the end 1652 -1 samtools index: "sample20bowtiemapping_sorted.bam" is corrupted or unsorted @Poshi

samtools faidx file.bam file.bai makes no sense. samtools faidx indexes fasta files, not bams or bais (bai are the indexes of bam files).

I'm not sure about the new error you are getting. It looks like the execution of some samtools index, but no samtools index is in the code. Can you update the question with the new information? Put the new code, tell us the line that fails and put also the error message.

Okay, Ive updated everything with the new error and new code  @Poshi

Which is the command that fails? Which one is giving that message? Again, I'm seeing several errors. faidx have too many parameters (what are you trying to achieve?), view will complain about not having a bam input, and sort have no input file given and the shell will complain about an orphan dollar sign. Four errors in for different lines. Better go step by step and check that one step works before going to the next.

Try upgrading to the most recent version of samtools, this bug may have been squashed already.

I think a lot of the questions in the comments could be solved with a Minimal, Complete, and Verifiable example: https://stackoverflow.com/help/mcve

By now this will be of little interest to the original poster, but for future reference:
$ samtools sort -n -T ~/[etc] -o sample${r}bowtiemapping_sorted.bam [etc]
$ samtools index sample${r}bowtiemapping_sorted.bam
[E::hts_idx_push] NO_COOR reads not in a single block at the end 1652 -1.

This error message indicates (admittedly somewhat obscurely) that the reads are not sorted by coordinate — in particular, there are mapped reads appearing after the unplaced unmapped reads (those with RNAME *) that should be at the very end of the file.
This is because samtools sort -n has been used to sort the reads by name instead. Remove -n to sort by position, which is what is needed to prepare a BAM file for indexing with samtools index.

Hi @JohnMarshall agreed but it certainly of interest to many of us.


text:
Fastest way to demultiplex bam file based on field
I have a bam file with aligned reads from multiple experiments. For each experiment I have a portion of the sequence identifier, thus for example

HWI-D00310:448:CCG0KANXX:7:1101:2570:1976,

and is different from other the other of single-run flow cells in the analysis.
I want to demultiplex the input bam file into different bam files depending on the experiment. My only concern is a fast methodology to demultiplex a bam using combinations of single-run, flow-cell, and lane. What I came up with (the fastest) is something like:
# Parallelized per chromosome. Assuming only 2 plates here.
printf '%s\n' $(seq 1 22) X Y MT \
| parallel 'samtools view input.bam {} \
  | tee >(grep -f patterns_experiment1.txt > experiment1_chr{}.sam.noheader) \
        >(grep -f patterns_experiment2.txt > experiment2_chr{}.sam.noheader)' \
&& seq 1 2 \
| parallel 'samtools view -H input.bam && cat experiment{}_*sam.noheader \
  | samtools view -1 - > input_experiment{}.bam'

Is there a better/faster way to do it?
If I understood correctly read groups should be specified prior to mapping? I started from the already mapped BAM unfortunately, but looking at the original FASTQs I noticed the pattern in which some parts of the sequence identifier were unique for one experiment vs the others.


text:
Codon usage analysis for whole genomes
I am new to bioinformatics. So if these questions seem you to a bit childish please forgive me.
I have two queries.

I am intending to perform a codon usage analysis followed by correspondence analyses for multiple microbial whole genomes of one bacterial species to find the association with the isolation source and ST type. Is this the right approach?
I have concatenated all the CDS in a single genome and joined all the genomes (average genes per genome is 4100). The problem is all the CDS when concatenated also have the stop codons in it. CodonW, famous for this analyses, cannot begin the analyses in the presence of stop codons. How can I remove stop codons?

Kindly suggest me a solution to the problem and concept.

This python script will remove all stop codons from your fasta file. I called it remove_stops.py
#!/usr/bin/env python3
import sys
from Bio import SeqIO
stop_codons = ["TAG", "TGA", "TAA", "UAG", "UGA", "UAA"]

fasta_it = SeqIO.parse(open(sys.argv[1]), 'fasta')
for fasta in fasta_it:
    name, sequence = fasta.description, str(fasta.seq)
    try:
        assert(len(sequence) % 3 == 0)
    except:
        raise IOError("length of sequence {} was not divisible by three".format(name))
    print(">{}".format(name))
    print("".join([sequence[i:i+3] for i in range(0, len(sequence) - 3 + 1, 3)
                   if sequence[i:i+3] not in stop_codons]))

You can run it like so:
python remove_stops.py stop_codons.fasta > filtered.fasta

This program won't work if there are any fasta records that do not have a length divisible by three. This indicates a sliding translation stop and these sequences should be trimmed individually. Also, all the genes must be in the 5'->3' direction.
It works by looking at each codon [sequence[i:i+3] for i in range(0, len(sequence) - 3 + 1, 3), then just makes sure they are not a stop codon: if sequence[i:i+3] not in stop_codons]


text:
levels of factors in the design have non-unique level names after make.names() is applied
When I try to import the data using DESeq2 package:
dds <- DESeqDataSetFromTximport(txi, sampleTable, ~Group + Cre)

I am getting the error:

Error in validObject(.Object) :
    invalid class “DESeqDataSet” object: levels of factors in the design have non-unique level names after make.names() is applied.
    best to only uobject letters and numbers for levels of factors in the design

Here is how my sampleTable looks like:
    ID      Group Cre
1  Bla_36   IL4   -
2  Bla_38   IL4   -
3  Bla_11   IL4   -
4  Bla_17   IL4   -
5  Bla_20   IL4   +
6  Bla_21   IL4   +
7   Bla_8   IL4   +
8  Bla_22   IL4   +
9   Bla_7   IL4   +
10 Bla_35   PBS   -
11 Bla_37   PBS   -
12 Bla_13   PBS   -
13  Bla_1   PBS   -
14  Bla_9   PBS   +
15 Bla_23   PBS   +
16 Bla_19   PBS   +
17 Bla_24   PBS   +

What am I doing wrong here? And how to make it working?

There's a good reason to not allow + and - in the columns. It's because once you do that you get coefficients that can then start with a + or -, which is illegal in R and most other languages. In R, legal variable names consist of:

A syntactically valid name consists of letters, numbers and the
       dot or underline characters and starts with a letter or the dot
       not followed by a number.

That doesn't really explain the error message you saw, of course. That can be solved by thinking about what R does to invalid character in names. Namely, it replaces them by a dot. So + and - get converted to X. and X., which are of course identical and are what's really leading to the error message your saw. If you had used +cre and - it would have worked.
As an aside, you'll want to think in terms of positive and negative Cre if you talk to a biologist, since that's the appropriate terminology for it.

Weirdly enough, + and - are apparently not allowed in the column. After mapping them to longer strings it started to work:
library(plyr)
cre_id = revalue(samples$cre, c('-' = 'minus', '+' = 'plus'))

sampleTable <- data.frame(
  ID = samples$id,
  Group = samples$group,
  Cre = cre_id
)

and it looks like that:
       ID   Group   Cre
1  Bla_36   IL4 minus
2  Bla_38   IL4 minus
3  Bla_11   IL4 minus
4  Bla_17   IL4 minus
5  Bla_20   IL4  plus
6  Bla_21   IL4  plus
7   Bla_8   IL4  plus
8  Bla_22   IL4  plus
9   Bla_7   IL4  plus
10 Bla_35   PBS minus
11 Bla_37   PBS minus
12 Bla_13   PBS minus
13  Bla_1   PBS minus
14  Bla_9   PBS  plus
15 Bla_23   PBS  plus
16 Bla_19   PBS  plus
17 Bla_24   PBS  plus

I do not understand though, why...


text:
Mapping RefSeq-ID or NCBI-ID (Entrez-ID) to EC number
I annotated my bacterial genomes using the new NCBI Prokaryotic Genome Annotation Pipeline and now, I want to annotate EC-numbers.
In the master annotation file (.gff), I get the RefSeq-ID of every successful annotation, e.g. RefSeq:YP_805528.1.
When I visit the respective NCBI page, I see that it sometimes is mapped to an EC number.
I'd like to automate this process: Convert every RefSeq-ID to an EC number. How can I do this?
(I may already have discovered the first step: I can map the RefSeq-ID to the NCBI-ID (Entrez-ID) using this NCBI conversion database: YP_805528.1 -> 4419225. But I don't see how to get from these databases to EC.)
Are there other ways? As far as I know, most people blast their genes to already-mapped databases like SwissProt and ExPASy ENZYME. But I don't like this approach because technically, the NCBI annotation pipeline already did this. Why compute the same thing twice?

You can use Entrez Direct for this as follows:
efetch -db protein -id YP_805528.1 -format gpc -mode xml | xtract -insd Protein EC_number
YP_805528.1     6.3.2.13


text:
What is the most compact data structure for canonical k-mers with the fastest lookup time?
edit: Results are current as of Dec 4, 2018 13:00 PST.
Background
K-mers have many uses in bioinformatics, and for this reason it would be useful to know the most RAM-efficient and fastest way to work with them programmatically. There have been questions covering what canonical k-mers are, how much RAM k-mer storage theoretically takes, but we have not yet looked at the best data structure to store and access k-mers and associated values with.
Question
What data structure in C++ simultaneously allows the most compact k-mer storage, a property, and the fastest lookup time? For this question I choose C++ for speed, ease-of-implementation, and access to lower-level language features if desired. Answers in other languages are acceptable, too.
Setup

For benchmarking:

I propose to use a standard fasta file for everyone to use. This program, generate-fasta.cpp, generates two million sequences ranging in length between 29 and 300, with a peak of sequences around length 60.
Let's use k=29 for the analysis, but ignore implementations that require knowledge of the k-mer size before implementation. Doing so will make the resulting data structure more amenable to downstream users who may need other sizes k.
Let's just store the most recent read that the k-mer appeared in as the property to retrieve during k-mer lookup. In most applications it is important to attach some value to each k-mer such as a taxon, its count in a dataset, et cetera.
If possible, use the string parser in the code below for consistency between answers.
The algorithm should use canonical k-mers. That is, a k-mer and its reverse complement are considered to be the same k-mer.

Here is generate-fasta.cpp. I used the command g++ generate_fasta.cpp -o generate_fasta to compile and the command ./generate_fasta > my.fasta to run it:
return 0;
//generate a fasta file to count k-mers
#include <iostream>
#include <random>

char gen_base(int q){
  if (q <= 30){
    return 'A';
  } else if ((q > 30) && (q <=60) ){
    return 'T';
  } else if ((q > 60) && (q <=80) ){
    return 'C';
  } else if (q > 80){
    return 'G';
  }
  return 'N';
}

int main() {
  unsigned seed = 1;
  std::default_random_engine generator (seed);
  std::poisson_distribution<int> poisson (59);
  std::geometric_distribution<int> geo (0.05);
  std::uniform_int_distribution<int> uniform (1,100);
  int printval;
  int i=0;
  while(i<2000000){
    if (i % 2 == 0){
      printval = poisson(generator);
    } else {
      printval = geo(generator) + 29;
    }
    if (printval >= 29){
      std::cout << '>' << i << '\n';
      //std::cout << printval << '\n';
      for (int j = 0; j < printval; j++){
        std::cout << gen_base(uniform(generator));
      }
      std::cout << '\n';
      i++;
    }
  }
  return 0;
}

Example
One naive implementation is to add both the observed k-mer and its reverse complement as separate k-mers. This is obviously not space efficient but should have fast lookup. This file is called make_struct_lookup.cpp. I used the following command to compile on my Apple laptop (OS X): clang++ -std=c++11 -stdlib=libc++ -Wno-c++98-compat make_struct_lookup.cpp -o msl.
#include <fstream>
#include <string>
#include <map>
#include <iostream>
#include <chrono>
//build the structure. measure how much RAM it consumes.
//then measure how long it takes to lookup in the data structure

#define k 29

std::string rc(std::string seq){
  std::string rc;
  for (int i = seq.length()-1; i>=0; i--){
    if (seq[i] == 'A'){
      rc.push_back('T');
    } else if (seq[i] == 'C'){
      rc.push_back('G');
    } else if (seq[i] == 'G'){
      rc.push_back('C');
    } else if (seq[i] == 'T'){
      rc.push_back('A');
    }
  }
  return rc;
}

int main(int argc, char* argv[]){
  using namespace std::chrono;
  //initialize the data structure
  std::string thisline;
  std::map<std::string, int> kmer_map;
  std::string header;
  std::string seq;
  //open the fasta file
  std::ifstream inFile;
  inFile.open(argv[1]);

  //construct the kmer-lookup structure
  int i = 0;
  high_resolution_clock::time_point t1 = high_resolution_clock::now();
  while (getline(inFile,thisline)){
    if (thisline[0] == '>'){
      header = thisline.substr(1,thisline.size());
      //std::cout << header << '\n';
    } else {
      seq = thisline;
      //now add the kmers
      for (int j=0; j< thisline.size() - k + 1; j++){
        kmer_map[seq.substr(j,j+k)] = stoi(header);
        kmer_map[rc(seq.substr(j,j+k))] = stoi(header);
      }
      i++;
    }
  }
  std::cout << "  -finished " << i << " seqs.\n";
  inFile.close();
  high_resolution_clock::time_point t2 = high_resolution_clock::now();
  duration<double> time_span = duration_cast<duration<double>>(t2 - t1);
  std::cout << time_span.count() << " seconds to load the array." << '\n';

  //now lookup the kmers
  inFile.open(argv[1]);
  t1 = high_resolution_clock::now();
  int lookup;
  while (getline(inFile,thisline)){
    if (thisline[0] != '>'){
      seq = thisline;
      //now lookup the kmers
      for (int j=0; j< thisline.size() - k + 1; j++){
        lookup = kmer_map[seq.substr(j,j+k)];
      }
    }
  }
  std::cout << "  - looked at " << i << " seqs.\n";
  inFile.close();
  t2 = high_resolution_clock::now();
  time_span = duration_cast<duration<double>>(t2 - t1);
  std::cout << time_span.count() << " seconds to lookup the kmers." << '\n';

}

Example output
I ran the above program with the following command to log peak RAM usage. The amount of time the lookup of all k-mers in two million sequences is reported by the program. /usr/bin/time -l ./msl my.fasta.
The output was:
 -finished 2000000 seqs.
562.864 seconds to load the array.
  - looked at 2000000 seqs.
368.734 seconds to lookup the k-mers.
     1046.94 real       942.38 user        78.96 sys
11680514048  maximum resident set size

So, the program used 11680514048 bytes = 11.68GB of RAM and it took 368.734 seconds to lookup the k-mers in two million fasta files.
Results
Below is a plot of the results from each user's answers.

try std::unordered_map instead of std::map. This should already give you quite a boost.

@PeterMenzel, thanks for your comment. What do you think about submitting this comment as an answer? I can implement it and add the results.

You can try std::unordered_map instead of std::map as a simple improvement (hopefully)

Can you be specific about which you care about more: compactness or speed? If both are equally important, then there will likely be no answer to this question.

Your generate-fasta.cpp does not compile. I also don’t understand its logic; why the multiple random distributions? What are you modelling? The gen_base function confuses me.

YMMV depending on your compiler. It worked for me with g++ Apple LLVM version 10.0.0 (clang-1000.10.44.2). The gen_base function uses a random int between 1-100 as input to select a base. The ranges in the function give it some AT bias, The Poisson distribution makes a bunch of sequences around 59 bp long, but the geometric distribution gives it a long tail of sequence lengths. In this case I was modeling a distribution of some sequences I need to look up kmers for.

You might mark the other answer as the correct one.

The question and the accepted answer are not about k-mer data structure at all, which I will explain in detail below. I will first answer the actual question OP intends to ask.
The simplest way to keep k-mers is to use an ordinary hash table. The performance is mostly determined by the hash table library. std::unordered_map in gcc/clang is one of the worst choices because for integer keys, it is very slow. Google dense, ska::bytell_hash_map and ska::flat_hash_map, tsl::robin_map and absl::flat_hash_map are much faster. There are a few libraries that focus on smaller footprint, such as google sparse and sparsepp, but those can be a few times slower.
In addition to the choice of hash table, how to construct the key is critical. For k<=32, the right choice is to encode a k-mer with a 64-bit integer, which will be vastly better than std::string. Memory alignment is also important. In C/C++, as long as there is one 8-byte member in a struct, the struct will be 8-byte aligned on x86_64 by default. Most C++ hash table libraries pack key and value in std::pair. If you use 64-bit keys and 32-bit values, std::pair will be 8-byte aligned and use 16 bytes, even though only 12 bytes are actually used – 25% of memory is wasted. In C, we can explicitly define a packed struct with __attribute__ ((__packed__)). In C++, probably you need to define special key types. A better way to get around memory alignment is to go down to the bit level. For read mapping, for example, we only use 15–23bp seeds. Then we have 18 (=64-23*2) bits left unused. We can use these 18 bits to count k-mers. Such bit-level management is quite common.
The above is just basic techniques. There are a few other tricks. For example, 1) instead of using one hash table, we can use 4096 (=2^12) hash tables. Then we can store 12 bits of k-mer information into the 4096 part. This gives us invaluable 12 bits in each bucket to store extra information. This strategy also simplifies parallel k-mer insertions as with a good hash function, it is rare to insert into two tables at the same time. 2) when most k-mers are unique, the faster way to count k-mers is to put k-mers in an array and then sort it. Sorting is more cache friendly and is faster than hash table lookups. The downside is that sort counting can be memory demanding when most k-mers are highly repetitive.

The other answer is spending considerable (probably the majority of) time on k-mer iteration, not on hash table operations. The program loops through each position on the sequence and then each k-mer position. For an $L$-long sequence, this is an $O(kL)$ algorithm. It has worse theoretical time complexity than hash table operations, which is $O(L)$. Although hash table operations are slow due to cache misses, a factor of k=29 is quite significant. Another issue is that all programs in the question and in the other answer are compiled without -O3. Adding this option brings the bytell_hash_map lookup time from 314s to 34s on my machine.
The C program at the end of my post shows the proper way to iterate k-mers. It is an $O(L)$ algorithm with a tiny constant. The program keeps track of both forward and reverse k-mers at the same time and update them with a few bit operations at each sequence position. This echoes my previous comment "You should not reverse complement the whole k-mer". On the same machine, the program looks up k-mers in 5.5s using 792MB RAM at the peak. This 6-fold (=34/5.5) speedup mostly comes from k-mer iteration, given that the hash table library in use is known to have comparable performance to bytell_hash_map.

#include <stdio.h>
#include <stdint.h>
#include "khash.h"

static inline uint64_t hash_64(uint64_t key)
{ // more sophisticated hash function to reduce collisions
  key = (~key + (key << 21)); // key = (key << 21) - key - 1;
  key = key ^ key >> 24;
  key = ((key + (key << 3)) + (key << 8)); // key * 265
  key = key ^ key >> 14;
  key = ((key + (key << 2)) + (key << 4)); // key * 21
  key = key ^ key >> 28;
  key = (key + (key << 31));
  return key;
}

KHASH_INIT(64, khint64_t, int, 1, hash_64, kh_int64_hash_equal)

unsigned char seq_nt4_table[128] = { // Table to change "ACGTN" to 01234
  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
  4, 0, 4, 1,  4, 4, 4, 2,  4, 4, 4, 4,  4, 4, 4, 4,
  4, 4, 4, 4,  3, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
  4, 0, 4, 1,  4, 4, 4, 2,  4, 4, 4, 4,  4, 4, 4, 4,
  4, 4, 4, 4,  3, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4
};

static uint64_t process_seq(khash_t(64) *h, int k, int len, char *seq, int is_ins)
{
  int i, l;
  uint64_t x[2], mask = (1ULL<<k*2) - 1, shift = (k - 1) * 2, tot = 0;
  for (i = l = 0, x[0] = x[1] = 0; i < len; ++i) {
    int c = (uint8_t)seq[i] < 128? seq_nt4_table[(uint8_t)seq[i]] : 4;
    if (c < 4) { // not an "N" base
      x[0] = (x[0] << 2 | c) & mask;                  // forward strand
      x[1] = x[1] >> 2 | (uint64_t)(3 - c) << shift;  // reverse strand
      if (++l >= k) { // we find a k-mer
        uint64_t y = x[0] < x[1]? x[0] : x[1];
        khint_t itr;
        if (is_ins) { // insert
          int absent;
          itr = kh_put(64, h, y, &absent);
          if (absent) kh_val(h, itr) = 0;
          tot += ++kh_val(h, itr);
        } else { // look up
          itr = kh_get(64, h, y);
          tot += itr == kh_end(h)? 0 : kh_val(h, k);
        }
      }
    } else l = 0, x[0] = x[1] = 0; // if there is an "N", restart
  }
  return tot;
}

#include <zlib.h>
#include <time.h>
#include <unistd.h>
#include "kseq.h"
KSEQ_INIT(gzFile, gzread)

int main(int argc, char *argv[])
{
  khash_t(64) *h;
  int i, k = 29;
  while ((i = getopt(argc, argv, "k:")) >= 0)
    if (i == 'k') k = atoi(optarg);
  h = kh_init(64);
  for (i = 1; i >= 0; --i) {
    uint64_t tot = 0;
    kseq_t *ks;
    gzFile fp;
    clock_t t;
    fp = gzopen(argv[optind], "r");
    ks = kseq_init(fp);
    t = clock();
    while (kseq_read(ks) >= 0)
      tot += process_seq(h, k, ks->seq.l, ks->seq.s, i);
    fprintf(stderr, "[%d] %.3f\n", i, (double)(clock() - t) / CLOCKS_PER_SEC);
    kseq_destroy(ks);
    gzclose(fp);
  }
  kh_destroy(64, h);
  return 0;
}

Thanks for updating your comments. Doesn't 'reduce collisions' in the hash function imply that collisions are possible? Collisions aren't OK for some applications. How often should this happen and under what circumstances?

@conchoecia you are misunderstanding what collisions mean in hash tables. Ordinary hash functions all have collisions. Perfect hash function doesn't, but constructing it is a non-trivial task and requires dedicated libraries. Perfect hash function is not always practical depending on applications.

What I think you mean by collisions -> "A collision means that two unique elements produce the same hash value. In this case for kmers, two unique DNA k-mer strings would produce the same hash value. So when counting 3-mers in the dataset AAAT, if AAT and AAA have the same hash value, our data structure would show that the k-mers AAT occurs twice and AAA occurs twice."

@conchoecia Hash table keeps the actual keys, not only their hashes. AAT and AAA are distinct keys and will only be counted once. Collision is something that every standard library has to resolve. Every implementation here has collisions. Please read the wiki page to understand how hash table works. PS: collisions have nothing to do with the results. They only affect performance.

The more I go over your code the more I am learning. This is amazing. Thank you so much.


text:
Problem with SPAdes installation - no bin folder
I have been trying to download SPAdes for an important project I wanted to do, however, I have encountered a small problem while doing so. It appears that there is no 'bin' folder in my SPAdes installation file.
I have used 7-Zip as my main method of unpacking the file (As I am using Windows, I had to download the actual source code rather than a nice executable), and it came in the .tar.gz format, which I believe 7-Zip is rather good with.
A few pictures detailing my problems are shown below.
I hope I can get a little help on this matter. It is slightly time-sensitive, so I would really appreciate a little bit of assistance shortly. Thank you in advance!

As seen above, the first image is the command I used, the usual test command, the second is the file structure I got (Notice the absence of the 'bin' file), and the third is the error I received from running the spades.py --test command.

I see that it supports Linux and MacOS. It appears you are using Windows. Is there a reason you believe that Windows is supported? See https://github.com/ablab/spades/issues/156

No, not at all. I just thought the source code would work on Windows, as there are options for Linux and MacOS downloadables.

Well, sure, 7-Zip will unzip a tar.gz file, but it won't install software for you!

Simply unpacking the source code won't actually install much, you'll need to follow its installation instructions. Note that this might be possible on Windows if you use new bash command line, but in general don't expect any bioinformatics tools to work on windows.


text:
Integrative analysis of omics studies using machine learning
I would like to use public omics datasets (ChIP-seq, RNA-seq, and ATAC-seq) from different studies to do an integrative analysis as follow:

Normalise samples, within each type of omics, from different public datasets.
Convert the normalised values into a uniform scale to make the comparison between ChIP-seq, RNA-seq and ATAC-seq possible.
Feed the normalised uniformed values into machine learning to infer one feature (e.g. RNA expression) from other features (e.g. TF or histone marks ChIP-seq).

How could I integrate these data-sets to predict one group of variables from another one?

This question is very general and opinion based.  Could you please make the question more specific

What is your goal? Do you have some data at hand? Would the data of the different omic come from the same sample or from different subjects? Would the data be quite homogeneous or not? How many samples do you have for each dataset ?

These are mouse samples, same background, so I am assuming homogeneity.
I have at least 2 replicates per samples, where sample A --is a precursor to --> B --> C.
My goals are:
a) to describe the dynamic of different histone marks during differentiation (A -> B -> C).
b) to predict gene expression from these histone marks.

Could you further elaborate that @Firas? How many samples do you have (6)? Do you have two replicates per A, B and C or more? Are they technical replicates or biological replicates? With a bit more of information I could possibly answer it

If you are looking for transcription regulation, then the following article will give you a basic idea. https://www.cell.com/developmental-cell/fulltext/S1534-5807(16)30002-8

@Llopis Many thanks, I have 5 histone chip-seq, atac-seq and RNA-seq in 2 technical replicates for each of A, B and C.

@arupgsh great reference, I will look into it. Many thanks.

The steps you describe are correct. For step 2 it is usually normalized to mean 0 and variance 1. However the "machine learning" part is important.
Having several samples being technical replicates will make the integration task easier. However, you have too few samples to make any good prediction. At most I would describe it as an exploratory analysis.
I recommend using regularized canonical correlation analysis for this kind of multi-omics analysis. The method seeks which are the features of one dataset that correlate (or covariate, depending on a parameter) with the variables of other datasets. Then you could use those features to predict the position of the other.
There are several implementations of this methods, the original in the RGCCA package, and mixOmics for other procedures with cross-validations also STATegRa has other related methods.


text:
How to convert Bed file to fasta file?
I have Bed file containing start and end of a sequence, and I need to convert it to fasta format, any recommendations?

Please [edit] your question and give us more details. Bed files define regions, do you have a fasta file of the target regions? Is this for a genome? Do you have that genome?

bedtools getfasta can do this. See documentation.

You could start with a folder of samtools faidx-indexed FASTA files for your genome of interest:
$ mkdir /genomes/hg38 && cd /genomes/hg38
$ wget ftp://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/*.fa.gz
$ for fn in `ls *.fa.gz`; do gunzip $fn; done
$ for fn in `ls *.fa`; do samtools faidx $fn; done

Then use that folder with the following script:
#!/usr/bin/env perl

use strict;
use warnings;
use Getopt::Long;

#
# bed2faidxsta.pl
# --
# Reads BED data from standard input, writes FASTA to standard output.
#
# Dependent on samtools and indexed FASTA data files located in $fastaDir
# variable. Set --fastaDir=dir to set custom directory containing a source
# of per-build, bgzip-compressed FASTA and associated index (fa.gz.fai)
# files, or leave unset to use data in current working directory. Use the
# --fastaIsUncompressed option if the FASTA files are not compressed.
#

# test if samtools is available
`samtools --version` || die "Error: The samtools application is required to run this script. Try 'module add samtools' or install a local copy of samtools.\n";

# default FASTA input is current working directory
my $fastaDir = `pwd`; chomp $fastaDir;
# default is to assume input coordinates use zero-based index scheme
my $oneBased;
# default is to leave IDs alone
my $useIDPrefixAsStrand;
# default is to assume FASTA files are bgzip-compressed
my $fastaIsUncompressed;

GetOptions ('fastaDir=s' => \$fastaDir, 'oneBased' => \$oneBased, 'useIDPrefixAsStrand' => \$useIDPrefixAsStrand, 'fastaIsUncompressed' => \$fastaIsUncompressed);

if (! -d $fastaDir) { die "Error: FASTA directory does not exist\n"; }

while (<STDIN>) {
    chomp;
    my ($chr, $start, $stop, $id, $score, $strand) = split("\t", $_);
if (!defined($chr) || !defined($start) || !defined($stop)) { die "Error: No chromosome name, start or stop position defined\n"; }
    if (!defined($id)) { $id = "."; }
if (!defined($score)) { $score = "."; }
if (!defined($strand)) { $strand = "+"; } else { $strand = substr($strand, 0, 1); }
# adjust coordinates to one-based index, if necessary
my ($queryChr, $queryStart, $queryStop) = ($chr, $start, $stop);
if (!$oneBased) {
        $queryStart++;
}
# adjust strand if required
if ($useIDPrefixAsStrand) {
        $strand = substr($id, 0, 1);
    }
    # lookup
    my $queryFn = "$fastaDir/$chr.fa.gz";
if ($fastaIsUncompressed) {
        $queryFn = "$fastaDir/$chr.fa";
}
my $queryKey = "$queryChr:$queryStart-$queryStop";
my $queryResult = `samtools faidx $queryFn $queryKey`; chomp $queryResult;
# linearize result
my @lines = split("\n", $queryResult);
    my @seqs = @lines[1..(scalar @lines - 1)];
    my $seq = join("", @seqs);
# handle reverse-stranded elements
if ($strand eq "-") {
        $seq = rc_sequence($seq);
    }
    # print to standard output
    my $header = ">".join(":",($chr, $start, $stop, $id, $score, $strand));
print STDOUT $header."\n".$seq."\n";
}

sub rc_sequence {
    my $seq = shift @_;
my $reverse_complement = reverse($seq);
$reverse_complement =~ tr/ACGTacgt/TGCAtgca/;
    return $reverse_complement;
}

Then:
$ perl bed2faidxsta.pl < in.bed > out.fa

This approach can be useful with stranded input, nonstandard genomes, or just to do things with a bit more knowledge of inputs.

Or, since you already have samtools, the poster can use samtools faidx with a file of regions.  It won't be hard to get a bedfile into the right format: http://www.htslib.org/doc/samtools.html

UCSC 'kent' command line tool:
faToTwoBit -bed=yourCoordinates.bed genomeSequence.2bit output.fa

Binary for this command and others available at: UCSC kent command line tools

I think you mean twoBitToFa. Thanks anyway for pointing me in that direction.


text:
List of human transcription factors
I'm looking for a list of human transcription factor (TF) genes.
Where can I find a TF database or how can I retrieve this list from other databases?

Lambert et al. (2018, Cell. 172(4):650-665) provide a manually curated list of human TFs at this website.


text:
Build a Pubmed query given long gene list
I would like to search the literature using the Pubmed search engine for a set of genes that are associated with a brain region of interest. Unfortunately, there are a lot of them. It's possible there is no solution if there is a character limit on the query. I think it would be easiest to do this with 'nix tools but I'm open to R libraries or Python implementations too.
The list is 908 genes long. It's from the Gensat annotation search engine, with the filters Structure IS olfactory bulb & Expression level IS moderate to strong signal & Expression pattern IS region-specific.
The list looks like, scraped from the dropdown list in the HTML source with a regular expression:
Aatf
Abcb1b
Abhd3
Abhd4
Abhd5
...

The Pubmed query should look like:
(((((((Aatf) OR Abcb1b) OR Abhd3) OR Abhd4) OR Abhd5) OR Abl2) AND olfactory bulb

but with 908 genes separated by OR and parentheses.
Bonus points if you find/write a utility which adds gene synonyms to the query!

Perhaps you could get help on academia.stackexchange.com better than here. But what do you expect from such a long query? Could it work for you to search  this geneset/pathway on a pathway database and find the associated literature?

Mainly to subset the papers which may be relevant to these genes in the context of this tissue, then to manually comb through them to identify cell types they might be associated with.

Then your f"final" objective is to find which cell lines are associated with this genes? [Edit] the question to add this it might be someone knows a database of cell lines or something alike

That is out of the scope of my question. The goal is manual verification of cell types from sources in the literature. I do appreciate that there may be alternative solutions with databases already annotated.

After getting it out on paper (so to speak) I was able to accomplish what I wished with bash:
#!/bin/bash

gene_list=($(cat ./markers_clean.txt))
query=${gene_list[@]:0:1}
parens="("
for gene in ${gene_list[@]:1}
    do
        query="${query} OR ${gene})"
        parens="${parens}("
    done
echo "${parens}${query} AND olfactory bulb)"

With the output being:
((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((Aatf OR Abcb1b) OR Abhd3) OR Abhd4) OR Abhd5) OR Abl2) OR Ablim2) OR Acaa2) OR Ackr3) OR Acsbg1) OR Acvr1c) OR Acvr2a) OR Adam11) OR Adam23) OR Adcy1) OR Adcyap1) OR Adgrg1) OR Adora2a) OR Adra2a) OR Adra2c) OR AF529169) OR Agrp) OR Agt) OR Agtr1a) OR Agtr2) OR Ak4) OR Akap12) OR Akap5) OR Akr1c18) OR Aldh1l1) OR Aldh7a1) OR Aldoc) OR Als2) OR Amfr) OR Amigo1) OR Amigo2) OR Angpt2) OR Ank) OR Ank2) OR Ankrd34b) OR Api5) OR Apod) OR Arc) OR Arf6) OR Arhgap33) OR Arhgdig) OR Arhgef6) OR Arhgef7) OR Arl6ip1) OR Arpp21) OR Arvcf) OR Arx) OR Ascl1) OR Aspm) OR Atg4c) OR Atic) OR Atoh1) OR Atp1a1) OR Atp1b1) OR Atp2a2) OR Atp5a1) OR Atp5j2) OR Atp6v1b2) OR Atpif1) OR Avp) OR Avpr1a) OR Avpr1b) OR B2m) OR B3gat1) OR Bahcc1) OR Baiap2) OR Bcar1) OR Bdnf) OR Becn1) OR Bend5) OR Bmp1) OR Bmp6) OR Bmp7) OR Bnip3) OR Bsx) OR Btbd17) OR Btg1) OR Btg2) OR Btg3) OR C1ql2) OR C1ql3) OR C2cd2l) OR C33) OR Cacna1h) OR Cacna1i) OR Cacnb1) OR Cacng3) OR Cacng5) OR Cadps) OR Calb1) OR Calb2) OR Calca) OR Calcb) OR Calcr) OR Camk2a) OR Camkk1) OR Camkk2) OR Camkv) OR Card6) OR Cartpt) OR Casp9) OR Casr) OR Cbln1) OR Cbln2) OR Cbln4) OR Cbx5) OR Ccdc163) OR Ccdc184) OR Ccdc187) OR Ccdc3) OR Cck) OR Cckbr) OR Ccna2) OR Ccnd1) OR Ccne1) OR Ccne2) OR Ccsap) OR Cd63) OR Cdc42bpg) OR Cdc42ep4) OR Cdc6) OR Cdh1) OR Cdh11) OR Cdh3) OR Cdh6) OR Cdh9) OR Cdhr1) OR Cdk17) OR Cdk4) OR Cdk5) OR Cdk5r1) OR Cdk5r2) OR Cdk6) OR Cdkn1a) OR Cdkn1b) OR Cdkn2a) OR Cdkn2c) OR Cdkn2d) OR Cdon) OR Cdr2) OR Celf6) OR Cemip) OR Chd7) OR Chn2) OR Chrm1) OR Chrm4) OR Chrna2) OR Chrna3) OR Chrna4) OR Chrna5) OR Chrna7) OR Chrnb2) OR Chrnb4) OR Chst15) OR Chst2) OR Cic) OR Cisd1) OR Cited1) OR Ckmt1) OR Clcn1) OR Cldn1) OR Clu) OR Cnr1) OR Cntfr) OR Cntn2) OR Cntnap1) OR Coch) OR Col6a1) OR Colgalt2) OR Cpe) OR Cpne6) OR Crabp1) OR Creb1) OR Crh) OR Crhr1) OR Crhr2) OR Crim1) OR Crip2) OR Crk) OR Crym) OR Csnk1g3) OR Ctgf) OR Ctxn1) OR Cx3cl1) OR Cxcl11) OR Cxcl12) OR Cxcr4) OR Cxxc4) OR Cyb561) OR Cyp19a1) OR Cyp4v3) OR Cyr61) OR D8Ertd738e) OR Dbi) OR Dbp) OR Dcdc2a) OR Dct) OR Dctn1) OR Dcx) OR Ddit4l) OR Dgcr14) OR Dgcr6) OR Dhcr24) OR Dhrs3) OR Dhx57) OR Dkk3) OR Dlg4) OR Dlk1) OR Dll1) OR Dlx1) OR Dlx2) OR Dlx5) OR Dlx6) OR Dmbx1) OR Dnajb1) OR Dnajc12) OR Dnmt1) OR Dnmt3a) OR Doc2g) OR Drd1) OR Drd2) OR Drd3) OR Dsp) OR Dtx4) OR Dtymk) OR Dusp1) OR Dusp14) OR Dusp16) OR Dusp3) OR Dusp5) OR Dyrk1a) OR Dyx1c1) OR Dzip1) OR E2f8) OR Ebf2) OR Ebf3) OR Eef1a1) OR Efna1) OR Efnb1) OR Efnb2) OR Efr3a) OR Efs) OR Egr1) OR Eif1ax) OR Eif3h) OR Eif4e3) OR Eif4g2) OR Elavl2) OR Elavl3) OR Elavl4) OR Elfn1) OR Ell3) OR Emp1) OR Emx1) OR En2) OR Enah) OR Enc1) OR Eno1) OR Eno2) OR Entpd4) OR Eomes) OR Epha4) OR Epha7) OR Epha8) OR Eps8l2) OR Erp29) OR Esr2) OR Esyt3) OR Evx1) OR F2r) OR Fam171b) OR Fam81a) OR Fam84a) OR Fbxl5) OR Fbxo32) OR Fezf1) OR Fezf2) OR Fgd5) OR Fgf1) OR Fgf1) OR Fgf13) OR Fgf15) OR Fgfr1) OR Fign) OR Fkbp2) OR Fkbp3) OR Flrt1) OR Flrt2) OR Fn1) OR Fos) OR Foxi1) OR Foxi2) OR Foxm1) OR Foxp4) OR Frat1) OR Fryl) OR Frzb) OR Fxyd6) OR Fxyd7) OR Fzd2) OR Fzd3) OR Fzd5) OR Fzd8) OR Gaa) OR Gabbr1) OR Gabrb2) OR Gabre) OR Gabrg1) OR Gabrg2) OR Gad2) OR Gal) OR Galnt1) OR Galnt5) OR Galnt9) OR Galr3) OR Gap43) OR Gas1) OR Gbx2) OR Gck) OR Gcm2) OR Gdf1) OR Gfap) OR Gja1) OR Gjd2) OR Glce) OR Gldc) OR Glrb) OR Gls) OR Glud1) OR Gng1) OR Gng13) OR Gng3) OR Gng4) OR Gng5) OR Gng7) OR Gpc4) OR Gpm6a) OR Gpr12) OR Gpr139) OR Gpr151) OR Gpr165) OR Gpr21) OR Gpr27) OR Gpr37) OR Gpr68) OR Gpr83) OR Gpr88) OR Gprin3) OR Gpsm1) OR Gria2) OR Grid2ip) OR Grin2c) OR Grm1) OR Grm2) OR Grm4) OR Grp) OR Gsx2) OR Gtf2i) OR Gucy1a3) OR Hap1) OR Hbegf) OR Hcfc1) OR Hcn2) OR Hcn4) OR Hdac11) OR Hdc) OR Herc3) OR Hes6) OR Hey1) OR Hey2) OR Hgs) OR Hmgn2) OR Hmgn3) OR Hnrnpab) OR Hpcal1) OR Hrasls) OR Hrh2) OR Hsd11b1) OR Hspa12a) OR Hspb8) OR Htr1a) OR Htr1b) OR Htr1d) OR Htr3a) OR Htr3b) OR Htr4) OR Htr5b) OR Htr6) OR Icam5) OR Ict1) OR Id1) OR Id4) OR Idh3a) OR Ier5) OR Igf1) OR Igfbp2) OR Igfbp4) OR Igfbp5) OR Igsf9) OR Inpp5j) OR Insig2) OR Iqsec3) OR Irs4) OR Isl2) OR Islr2) OR Itgb3) OR Itm2b) OR Jag1) OR Jag2) OR Jak2) OR Jam3) OR Jph4) OR Jup) OR Kcnc1) OR Kcnc4) OR Kcnip2) OR Kcnj5) OR Kcnk1) OR Kcnq2) OR Kctd13) OR Kctd4) OR Kif11) OR Kif26a) OR Kit) OR Kitl) OR Klf1) OR Klf13) OR Kremen1) OR Krt19) OR Lamb1) OR Lamp5) OR Layn) OR Lbh) OR Ldb1) OR Lef1) OR Lfng) OR Lgals3) OR Lgi1) OR Lhx1) OR Lhx2) OR Lhx4) OR Lhx5) OR Lhx6) OR Lhx8) OR Limk1) OR Lix1) OR Lmo1) OR Lmo4) OR Lpar1) OR Lpl) OR Lrfn5) OR Lrig1) OR Lrig2) OR Lrp1) OR Lrp12) OR Lrp2) OR Lrrk2) OR Lrrn1) OR Lrrn3) OR Lrrtm2) OR Lss) OR Ly6h) OR Lynx1) OR Lypd1) OR Lypd2) OR Lypd6) OR Lzts1) OR Mafb) OR MAGEA8) OR Man1a2) OR Map1b) OR Map1lc3a) OR Map2k6) OR Map4k5) OR Mapk1) OR Mapk8) OR Mapt) OR March4) OR Matk) OR Mb21d2) OR Mbd2) OR Mbd4) OR Mc3r) OR Mc4r) OR Mdk) OR Med12) OR Mef2c) OR Meis1) OR Met) OR Mfap2) OR Mfap4) OR Mfn1) OR Mmp17) OR Mmp2) OR Mpc1) OR Mst1r) OR Msx1) OR Mtap) OR Mturn) OR Mxi1) OR Mycn) OR Nae1) OR Napb) OR Nars) OR Ncor1) OR Ndn) OR Ndnf) OR Ndrg1) OR Ndufv1) OR Ndufv2) OR Nectin3) OR Nefh) OR Nefl) OR Nefm) OR Neurod1) OR Neurod4) OR Neurod6) OR Neurog1) OR Neurog2) OR Nfil3) OR Nhlh1) OR Nhp2) OR NKAIN2) OR Nkain3) OR Nkx2) OR Nlgn3) OR Nmb) OR Nmbr) OR Nmu) OR Nnat) OR Nog) OR Nomo1) OR Nos1) OR Notch1) OR Notch2) OR Nov) OR Nova1) OR Npc1) OR Npc2) OR Nptx1) OR Nptx2) OR Npy) OR Npy1r) OR Npy2r) OR Npy5r) OR Nr2f2) OR Nr4a1) OR Nrep) OR Nrgn) OR Nrip3) OR Nrp1) OR Nrtn) OR Nrxn2) OR Nsf) OR Nsg1) OR Nsmf) OR Ntn1) OR Ntng2) OR Ntrk1) OR Ntrk2) OR Ntrk3) OR Nts) OR Nudt2) OR Numb) OR Numbl) OR Odc1) OR Olfm1) OR Olfr734) OR Olig2) OR Oma1) OR Ophn1) OR Opn3) OR Opn4) OR Oprk1) OR Otof) OR Otx2) OR Oxt) OR Oxtr) OR P2rx7) OR Pacsin3) OR Pak1) OR Pard6a) OR Pcbd1) OR Pcdh2) OR Pcdh9) OR Pcnt) OR Pcp4) OR Pcsk9) OR Pde1b) OR Pde1c) OR Pdp1) OR Pdyn) OR Peg3) OR Penk) OR Per2) OR Pex5l) OR Pgrmc1) OR Phactr1) OR Phactr4) OR Phox2b) OR Pias3) OR Pik3r1) OR Pirt) OR Pkd1) OR Plagl1) OR Plcxd2) OR Pld3) OR Plekhg1) OR Plk2) OR Plk3) OR Plp1) OR Plxdc1) OR Plxna1) OR Plxnc1) OR Pmp22) OR Pnmt) OR Pnoc) OR Podxl2) OR Polr1d) OR Polr3e) OR Pomc) OR Pou2f1) OR Pou3f4) OR Ppfia4) OR Ppm1e) OR Ppp1ca) OR Ppp1r14c) OR Ppp1r16b) OR Ppp1r17) OR Ppp1r1a) OR Ppp1r1b) OR Ppp1r2) OR Ppp1r7) OR Ppp2r2a) OR Ppp4c) OR Ppp6c) OR Prkacb) OR Prkcd) OR Prkcq) OR Prmt2) OR Prnp) OR Prodh) OR Prok2) OR Prokr2) OR Prox1) OR Prss12) OR Prtg) OR Psen1) OR Psen2) OR Ptch2) OR Pten) OR Ptf1a) OR Ptgds) OR Ptn) OR Ptpn2) OR Ptpn5) OR Ptprf) OR Pus1) OR Pygb) OR Rab37) OR Rac1) OR Ramp3) OR Rap1a) OR Rasa3) OR Rasal1) OR Rasgrp1) OR Rasgrp2) OR Rasl11b) OR Rbp1) OR Rcan2) OR Rcn1) OR Rel) OR Rest) OR Rfng) OR Rgs12) OR Rgs16) OR Rgs2) OR Rgs3) OR Rgs4) OR Rgs5) OR Rgs8) OR Rgs9) OR Rhob) OR Rhou) OR Rif1) OR Rlbp1) OR Rnd1) OR Rnd2) OR Rnd3) OR Rorb) OR Rpl18) OR Rpl31) OR Rprm) OR Rraga) OR Rtn4rl1) OR Rubcn) OR Rufy3) OR Runx1) OR Rxfp3) OR S1) OR Satb2) OR Scg2) OR Scg3) OR Scn3b) OR Scrn1) OR Scrt1) OR Scube1) OR Scube2) OR Sdc3) OR Selenom) OR Selenow) OR Sema3a) OR Sema3b) OR Sema3c) OR Sema3e) OR Sema3f) OR Sema4a) OR Sema6d) OR Sema7a) OR Sept3) OR Sept6) OR Serf1) OR Serpina9) OR Serpini1) OR Sertm1) OR Sez6) OR Sez6l2) OR Sf1) OR Sfrp1) OR Sfrp2) OR Sfrp5) OR Sim1) OR Sirpa) OR Sirt2) OR Skil) OR Sla) OR Slc1) OR Slc12a5) OR Slc13a3) OR Slc17a6) OR Slc17a7) OR Slc18a2) OR Slc1a1) OR Slc1a2) OR Slc1a3) OR Slc1a6) OR Slc22a3) OR Slc25a12) OR Slc32a1) OR Slc33a1) OR Slc41a3) OR Slc6a1) OR Slc6a12) OR Slc6a2) OR Slc6a3) OR Slc6a6) OR Slc6a9) OR Slc7a14) OR Slc8a2) OR Slc9a3r1) OR Slit1) OR Slitrk2) OR Slitrk3) OR Slitrk4) OR Slitrk5) OR Slitrk6) OR Smad6) OR Smarce1) OR Smoc2) OR Snap25) OR Snap91) OR Sncg) OR Snhg11) OR Snx17) OR Sort1) OR Sox1) OR Sox11) OR Sox2) OR Sox4) OR Sox8) OR Sox9) OR Sp4) OR Sprn) OR Spry1) OR Sptan1) OR Sqle) OR Srgap3) OR Srsf7) OR Sstr1) OR Sstr2) OR St8sia5) OR Stard5) OR Stau1) OR Strip2) OR Strn) OR Stum) OR Stx1a) OR Sufu) OR Sulf2) OR Sun2) OR Suv39h1) OR Sv2a) OR Sv2b) OR Syndig1l) OR Syngr3) OR Synj2) OR Syt13) OR Syt17) OR Syt2) OR Syt3) OR Syt4) OR Syt6) OR Syt7) OR Tac1) OR Tac2) OR Tacr3) OR Taf13) OR Tagln3) OR Tbc1d8) OR Tcf3) OR Tgif1) OR Thbs1) OR Thbs2) OR Thbs3) OR Thbs4) OR Them6) OR Thra) OR Thrsp) OR Tigar) OR Tle1) OR Tle3) OR Tmeff1) OR Tmem35a) OR Tmem59l) OR Tnc) OR Tor1b) OR Tpbg) OR Tpd52l1) OR Tph2) OR Tppp3) OR Traf4) OR Trappc3) OR Trf) OR Trh) OR Trhr) OR Trib2) OR Tsc22d2) OR Tspan14) OR Tspan7) OR Tuba1b) OR Tubb3) OR Tubb4a) OR Ube3a) OR Ubxn11) OR Uchl1) OR Ucp2) OR Uhmk1) OR Uncx) OR Uso1) OR Vat1) OR Vcan) OR Vdac2) OR Vegfa) OR Vim) OR Vip) OR Vipr2) OR Vldlr) OR Vmn1r2) OR Vsnl1) OR Vsx2) OR Washc2) OR Wdr46) OR Wfs1) OR Wif1) OR Wnt2b) OR Wnt5b) OR Wnt7a) OR Wnt8b) OR Wrb) OR Wscd1) OR Wwc1) OR Ypel4) OR Zbtb18) OR Zbtb2) OR Zcchc12) OR Zfp5) OR Zfp532) OR Zic1) OR Zic2) OR Zic3) OR Zic4) OR Zic5) OR Zwint) AND olfactory bulb)

As far as I can tell, this does not surpass any character limits on Pubmed.


text:
Collapse cell barcodes distribution within 1 Hamming distance
I have a barcode distribution from single-cell data, e.g:
11612552 TCCTGAGCACTGCATAACTCAA
9349711 GCTACGCTACTGCATAAGTCCA
8343678 CAGAGAGGCTAAGCCTGCACAT
8161950 CGTACTAGTCTCTCCGCGGCTA
8102383 TCCTGAGCGTAAGGAGCAGATC
7110298 AAGAGGCAGTAAGGAGATAATC
6626630 TAAGGCGAACTGCATAACATGT
6390489 CGAGGCTGTATCCTCTACTAGA
6210446 AGGCAGAATCTCTCCGAGTGCT
5985219 TAAGGCGATCTCTCCGTACGTG
...

where the first column is the number of reads associated with the cell barcode in the second column.
I would like to collapse each set of cell barcodes within 1 Hamming distance (i.e. allowing one mismatch between barcodes), keeping only the barcode with the highest read count and adding the read counts of the duplicates to it.
I was thinking to use awk and grep to filter the barcode list in an iterative procedure using regular expressions, but I think this would require too much time. I can share my code if needed.
Is there a way to do it using unix tools or already published pipelines (e.g. UMI tools)?
This is my code at the moment (it is not optimized though):
#!/bin/awk -f

{
    matching = 0;
    for (r in regexs) {
        if ($2 ~ regexs[r]) {
        matching = 1;
        counts[r] += int($1);
        break;
    }
}
if (!matching) {
    regex = "";
    for (i = 1; i <= length($2); i++) {
        l_regex = substr($2, 1, i-1) "." substr($2, i+1);
        if (length(regex) == 0)
            regex = l_regex;
        else
            regex = regex "|" l_regex;
    }
    regexs[$2] = regex;
 counts[$2] = int($1);
    }
    # so we know that 100 records have been processed
    if (!(NR % 100)) printf "." > "/dev/stderr";
}
END {
    for (c in counts) {
        print counts[c], c | "sort -k1 -rn";
    }
}

Doesn't umitools do this for you?

The problem with UMI tools (as I understood from here) is that it runs (whitelist) with fastqs, while I have already the count distribution. As of now, I haven't found a way to use directly the distribution with UMI tools.

So just go back to the fastq files and make your life easier.

It would be awesome, but before generating the count distribution, the analysis pipeline I'm using filters out reads based on multiple and supplemental alignments. Additionally, due to the size of the data I had to run everything through pipes, and the first step is the concatenation of the original fastqs, which were provided by the sequencing center already demultiplexed on a portion of cell barcode. My objective for this question is to collapse the distribution after filtering and mapping, for which I proposed a sample code, in the future I'll surely use UMI tools or similar from the start.

Also, if this is 10x, or Drop-seq, check out alevin, which is an all in one pipeline that includes dealing with multiple and supplemental alignments in an all-in-one pipeline that goes from fastq to counts matrix.

This doesn't use standard unix tools, but rather the API for UMI tools:
#!/usr/bin/env python
from umi_tools import network

cb = []  # cell barcodes
cbc = dict()  # cell barcode counts
for line in open("UMIs.txt"):
    cols = line.strip().split("\t")
    cols[1] = bytes(cols[1], "ascii")
    cb.append(cols[1])
    cbc[cols[1]] = int(cols[0])

uc = network.UMIClusterer()
CBclusters = uc(cb, cbc, 1)
cbFinal = dict()
for l in CBclusters:  # This is a list with the first element the dominant barcode
    cbFinal[l[0]] = 0
    for x in l:  # Iterate over all barcodes in a cluster
        cbFinal[l[0]] += cbc[x]
print(cbFinal)

If you give it input in UMIs.txt of the following form (note the last two lines differ from the first with a distance of 1):
11612552    TCCTGAGCACTGCATAACTCAA
9349711 GCTACGCTACTGCATAAGTCCA
8343678 CAGAGAGGCTAAGCCTGCACAT
8161950 CGTACTAGTCTCTCCGCGGCTA
8102383 TCCTGAGCGTAAGGAGCAGATC
7110298 AAGAGGCAGTAAGGAGATAATC
6626630 TAAGGCGAACTGCATAACATGT
6390489 CGAGGCTGTATCCTCTACTAGA
6210446 AGGCAGAATCTCTCCGAGTGCT
5985219 TAAGGCGATCTCTCCGTACGTG
20  TCCAGAGCACTGCATAACTCAA
2   TCCTGAGCACTGCATAACTCNA

Then it will produce:
{b'TAAGGCGATCTCTCCGTACGTG': 5985219,
 b'TAAGGCGAACTGCATAACATGT': 6626630,
 b'GCTACGCTACTGCATAAGTCCA': 9349711,
 b'CGTACTAGTCTCTCCGCGGCTA': 8161950,
 b'TCCTGAGCGTAAGGAGCAGATC': 8102383,
 b'CAGAGAGGCTAAGCCTGCACAT': 8343678,
 b'TCCTGAGCACTGCATAACTCAA': 11612574,
 b'AGGCAGAATCTCTCCGAGTGCT': 6210446,
 b'AAGAGGCAGTAAGGAGATAATC': 7110298,
 b'CGAGGCTGTATCCTCTACTAGA': 6390489}

Awesome, thanks! I didn't know you could use the API in this way.

Yeah, it doesn’t seem to be documented anywhere, but it’s python and the documentation in the source code is ok.

Sorry, we should get round to documenting this. Nice to see someone using it.

Note however, that this is using the "directional" algorithm. That is, the collapse will only happen if the two barcodes are within a hamming distance of 1 AND the counts for the second barcode are less that half those for the first. This was designed to work with UMIs rather than cell barcodes  There are other ways of deciding which collapses to make - see the paper for details. You can pass the method to the UMIClusterer constructor.

The difficulty with all these things comes with what to do about barcodes that are 1 edit distance away from more than one other barcodes in ways that form complex networks.
Cell barcodes are a bit different from UMIs. Partly this is because we expect them to be sparser. Partly because the consequences of inncorect error correction at the CB level are much worst than at the UMI level. The way that UMI-tools does its Cell Barcode error correction is:

Determine a whitelist of barcodes using the "knee" method
Iterate through the list of non-whitelist barcodes and check if they are within one edit distance of one and only one whitelist barcode. (this is done with umi_tools.umi_methods.getErrorCorrectMapping if you are interested).
Things that pass the above are collapsed into the whilelist barcode, else they are discarded.

The function umi_tools.umi_methods.getCellWhitelist takes a list of barcodes and their frequencies and returns a whitelist and (optionally) a mapping of which barcodes should be corrected to which others.


text:
Replace lowercase characters with -
I have an output from vcfutils.pl vcf2fq with specified minimal depth, and it means that nucleotides with not enough depth are lowercase.
I would like to change them to gaps. I could do it in higher language but I would like to know the solution in bash.
I have found simple solution to replace all lowercase to uppercase:
awk 'BEGIN{FS=" "}{if(!/>/){print toupper($0)}else{print $1}}' in.fna > out.fna

but I would like to replace lowercase with -. The first row need to be omitted. Could you help me? Thanks

The following sed command will do the trick.
sed -e '/^[^>]/ s/[a-z]/\-/g' in.fna > out.fna
Lines starting with > will be escaped and the lowercase characters will be replaced with - .

Thaks. I did try to do it without making new copy, but this works.

you can pipe this to vcfutils.pl vcf2fq output.

You could try it with sed.
Replace all lower-case chars with - in lines not starting with >:
sed -e '/^>/!s/[a-z]/-/g' in.fa

You can also use this solution where tr command from linux can solve your problem:
tail -n +2  in.fna | tr '[:lower:]' '-'