text:
Correlate DEGs from DESeq2, EdgeR and Limma results
I have a lists of DEGs identified by DESeq2, EdgeR and Limma.
I would like to correlate the the gene rankings in the lists to decide on a package to use in downstream analysis.
I am havig a few problems with this, would anyone be able to help?
Thanks

Hello Lyd, could you be a bit more specific? I think all these software have their advantages and disadvantages, it's impossible to give you a good advice without more information.
Also, did you tried to figure it by your own? What are the results of your research and why you are not convinced yet?

Im doing differential expresion analysis in R. Im comparing KD samples rto control samples using all 3 packages above. It was suggeested to me to do a correlation test between the results but Im no t sure what way to go about it. Not sure what the inputs  should be!

I would use the p-value/FDR which each method returns to rank the gene list for that method in order from 'most likely to be DE' to least likely.  You would then have three ranked lists of genes - you would put the actual rankings in for each gene - 1 for the top ranked, 2 for the next, etc.  You could then carry out, say Spearman's rank correlation on the lists of rankings to assess agreement between the methods.
You could also decide to keep all 3 methods, and accept genes as DE if, for example, at least 2 out of 3 methods give an acceptable p-value/FDR.
From my experience, correlation should be reasonably high at the top of the list, where the methods tend to agree more on the rankings, and when the correlation gets lower, further down the list, would be where I would start to question whether the genes are really differentially expressed.

This worked! Thank you :)


text:
How are residuals used as new phenotype?
I'm a novice trying to wrap my head around with fancy jargons. If I search residual, I understand that it represents the difference between the observed value of the dependent variable (y) and the predicted value (ŷ).
I'm unable to follow then, how are these residuals used as new phenotype?
I found similar question in the forum:  Adjusting phenotypes by regressing out covariates
The origin of my question is predixcan pipeline, where residuals are used as the new gene expression values after regressing input values on covariates.

The residual is the level of error in a regression model, the lower the residual the better the model. Residuals cannot equate to a phenotype, it is actually the opposite if your regression model is trying to investigate a given genotype/phenotype(s) model. The residual in this case is the amount of infomation which the genotype/phentype does not account for, i.e. there are other factors resulting in your observations which you have not captured in your model.
What I assume is that the regression model you are referring  is  investigating whether the residual is lowered by examining different gene combinations. The objective of a regression model is to lower the residual and I suspect that it is comparing the new residual against the old one. Therefore it is not being used as an input into a regression model but an assessment of its robustness.

Thank you very much for explaining this. Really appreciate it. :)

source: Stack Exchange

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

sort: votes


text:
High percentage of poly A sequences in 10X chromium R2 read
I'm currently analyzing two samples of eosinophil cells isolated from mouse lung and the samples are of very different quality.
According to the Cell Ranger summary 56% of the reads can be mapped to the transcriptome in the first sample and only 32% in the second sample. Also the number of cells identified in Cell Ranger are very different (5000 for sample 1 and 600 for sample 2).
When I look at the raw data in FastQC I see that the second sample has a lot of poly-A sequences in the R2 read that obviously can't be mapped to the transcriptome.
Sample 1:

Sample 2:

Also in many cases there is the following sequence in front of the poly-A:
AAGCAGTGGTATCAACGCAGAGTACATGGG
To me it looks as if most of the reads do not contain an insert and are primer dimers. Is this connected to the lower cell count in the second sample?
Maybe some of you have already seen something like this and have a suggestion on how to avoid such a problem in the future.
Thanks a lot

Hi @PPK, could you please add your solution as an answer to this question and "accept" it? That will help others who encounter the same problem, and stop the StackExchange Community bot from bumping this question from time to time to fish for more answers.

I actually figured out what was going on some time ago:
The over-represented sequence is the template switching oligo (TSO) that is used for the synthesis of the second strand (the example is for the v2 chemistry but the principle is the same for v3 and later).
This is a simplified illustration from the 10x manual:

As you can see the cDNA insert should be between the polyA capture probe and the TSO. In the second sample there are many cases with no insert and the TSO is followed by the poly A sequence. Actually the TSO should not be present in the final read after addition of Illumina adapters. But without any insert the fragmentation does not work and many reads still have this technical sequence.
So @merv was correct. Many reads contain no insert.
Why did this happen?
Why were there many reads without insert in this experiment and why are some samples more heavily affected than others? As I mentioned, this were eosinophil granulozytes and these fellas contain many proteases, DNAses and RNAses to destroy invading pathogens. Unfortunately, they also destroy mRNA in single cell experiments.
Interestingly, the eosinophily in the second sample received an activating treatment and probably contained even more or more active RNAses.

Answer from @merv, converted from comment:
Perhaps it's an insert size issue. Can you get the post library construction QC info? You might also see this in the BAMs, with the first showing a full bell-shaped peaks, and the latter truncating the peaks at the cleavage site. (Most 10X datasets I've seen have very few reads traversing the 3' end cleavage site.)

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


text:
Appropriate tool or algorithm for sloppy alignment of degenerate bases
I have an optimization problem where I have a degenerate nucleotide sequence I want to align against subsets of a reference genome (exons, specifically, to make the problem more tractable).
The sequence itself contains subregions of variable length.
The sequence looks like this, containing degenerate bases:
           L           C           R
5' - SSS......SSS|NNN.....NNN|WWW......WWW - 3'

Per IUPAC notation, S denotes C or G, N denotes any base, and W denotes A or T.
I am marking each block L, C, and R to broadly denote "left", "center", and "right" regions, oriented from 5' to 3' ("forward"-strand oriented).
The widths of the L and R blocks can be from 20 to 25 bases long. The width of the C block can be 14 to 16 bases long.
One programmatic way I could think to do this would be to brute-force a solution, generating all combinations of sequences of varied lengths, blastn-aligning each of such combinations to the sequence for the exon, and then looking for the best alignments among such sequences.
The output would be the scores for regions from sequences that align, to be able to say, from point X to point Y of the genome, a sequence here aligns better than other candidate sequences.
To avoid reinventing tools, and avoid brute-force methods, I'd like to know if there is a tool or even a wrapper around tools, which does this, in case it is already written, or another approach that is more elegant/efficient.

I believe mafft can align degenerate sequences.  From the documentation for version 7 at: https://mafft.cbrc.jp/alignment/software/anysymbol.html
Acceptable nucleotide symbols in the default mode:
a, c, g, t, u
Ambiguous nucleotides (r, y, w, s, k, m, d, v, h, b; IUPAC-IUB codes) can be used and are scored as:
score(r,a) = ( score(a,a) + score(g,a) ) / 2
score(r,g) = ( score(a,g) + score(g,g) ) / 2
score(r,t) = ( score(a,t) + score(g,t) ) / 2
score(r,c) = ( score(a,c) + score(g,c) ) / 2
...
score(r,y) = ( score(r,t) + score(r,c) ) / 2
...
score(r,r) = ( score(a,a) + score(g,g) ) / 2
score(y,y) = ( score(c,c) + score(t,t) ) / 2

MAFFT is designed for multiple sequence alignment, so will attempt to find a local alignment of homologous segments (using FFT), and then attempt to make a global alignment between the input sequences.

Can you expand your answer to include some support material?

I thought that a nucleotide blast could do the same. Just wondering if these tools also support querying against variable-width subsequences, or if wrappers around such tools can do this more efficiently in generalized cases. I have worked on parallelized analysis of combinations of substrings, in the meantime. It is brute-force but the scope of the problem seems tractable.

I'm not sure I quite understand your point about variable-width subsequences - as far as I understand, these tools tend to find the longest alignable subsequences they can from the two input sequences, so I guess if you want to test lost of different widths a wrapper would be the way to go.

@AlexReynolds I don't see why this isn't a 'normal' alignment problem where you could use blast or similar but with relaxed parameters (maybe permissive gap penalties to handle "variable-width" subsequences etc.). Perhaps a little more biological context would help us understand better too


text:
Creating a count table using R
I’m trying to develop a table from a series of lists generated using the Immunarch library to process TCR sequencing results. Each list is comprised of CDR3.aa (clone) information which are character strings and their count in a particular sample. The clones are short but vary between 7 and 20 (or more).  Each list has a header that identifies the sample sequenced. I have a list of 66 samples. Each sample can have several thousand clone strings. Not every clone is contained in every sample so the number of clones listed in the samples varies.  Here’s an example of the structure of a single sample list.
$PreAg_18_2
CDR3.aa           Count
CASSYGTAYTGELFF   1623
CASSRGDSDNSPLHF   1440
CASSREKAFF        1161
CSGMGALAKNIQYF     949
CSAYTGLSYEQYF      813
CASSLSLAVNSPLHF    634
CAIRDTPGSPQHF      574
CATGQVNTEAFF       555
CASSLKGQGGSPLHF    499
CASSYSRSPQPQHF     478

I want to combine the results in a single table showing the clone counts with all the clone strings listed on the y-axis and the sample ID listed on the x-axis. For example:
                10_pep_10_1     preAg_10_2      Dec_2_18_1  …...
CASSYGTAYTGELFF    1623         234             0
CASSRGDSDNSPLHF    1440         522             28
CASSREKAFF         1161         445             50
CSGMGALAKNIQYF      949         24              0
CASSYSRSPQPQHF      478         0               398
.
.

Currently, I'm working with R because that's what Immunarch is in, but I wonder if a) There's a better way to do this in Python, or b) If it's doable in r, then how do I generate the matrix matches.  The main problem is that not all the samples have the same number of clones, so I first have to separate each sample data before I can convert it to a data frame, and then how can I match the actual count with the actual clone.  I can't just merge the data frames.  Any suggestion is greatly appreciated.
Here's what i have for code right now.
Code:
    library(immunarch)
    library(stringr)
    library(plyr)
immdata = repLoad("/mnt/data/Development/Analysis_Scripts/input_files/")

all <- immdata$data

# Get list headers (names) and convert to df
sample.id <- names(all)
sample.id <- data.frame(sample.id)

# Get list of clones and filter for unique clones per list.
for (i in 1:length(all)){
    all[[i]]$Sample.ID<-names(all)[i]
    all[[i]]<-all[[i]][,c("CDR3.aa", "Clones")]
}

# all is the variable that contains the list of samples and clones.
all.split <- split(all, )

# make vector of all clones
all.clones <- unlist(all, use.names=FALSE)

# Removes clone repeats
all.clones.u <- unique(all.clones)

# convert list of clones to data frame
all.clones.u <- data.frame(all.clones.u)

At this point I have both the list of sample ids and the list of clones.  My problem is now matching each count with respective sample and  clone in the table.  Any suggestions?

It's hard to come up with a working code snippet if the workspace can't be recreated. Could you provide a minimal mock example of all that we can use? But based on what I think you want to do, I'd probably use a combination of left_join of the dplyr package combined with Reduce.

Thanks, Sebastian, I apologize for not providing a more clear picture of my current data, but I think you captured it correctly.  I'll give it a try.

You could give a try with dplyr's full_join() and purrr's reduce(). full_join() will make sure that you won't lose any data and reduce() will merge tables in your list from left to right, left being the first table of the list in the first iteration but would become merged table1+table2 in the second iteration and so on. This left to right "reduction" will also be helpful when renaming the columns of the resulting table.
Without your workspace I have attempted to replicate your data: There is a list of 3 data frames, each with 2 columns:
library(dplyr)
library(purrr)

df1 <- data.frame(CDR3.aa = c("A", "B", "C"),
                  Count = c(10,20,100))

df2 <- data.frame(CDR3.aa = c("A", "B", "D"),
                  Count = c(100,200,300))

df3 <- data.frame(CDR3.aa = c("D", "E", "F"),
                  Count = c(10,20,100))

my_list <- list(df1 = df1, df2 = df2, df3 = df3)

> my_list
$df1
  CDR3.aa Count
1       A    10
2       B    20
3       C   100

$df2
  CDR3.aa Count
1       A   100
2       B   200
3       D   300

$df3
  CDR3.aa Count
1       D    10
2       E    20
3       F   100

# "by" argument specifies which column to use when merging
my_table <- reduce(my_list, full_join, by = "CDR3.aa")

names(my_table) <- c("CDR3.aa", names(my_list))

> my_table
  CDR3.aa df1 df2 df3
1       A  10 100  NA
2       B  20 200  NA
3       C 100  NA  NA
4       D  NA 300  10
5       E  NA  NA  20
6       F  NA  NA 100

Clonotypes that are not present in table will be "NA".

Sounds like it could work.  I'll give it a try and report back.  Thanks.

haci, it worked very well.  I just had to convert the initial samples to data.frames so I could run reduce on it.  It's a very clean elegant solution.  Thank you very much for your help.

That can be easily done with MiXCR software.
If you first analyze the raw data using mixcr, which can be easily done with mixcr analyze command depending on your data structure:
e.g.
mixcr analyze generic-tcr-amplicon \
--species mmu \
--rna \
--rigid-left-alignment-boundary \
--floating-right-alignment-boundary C \
  input_R1.fastq.gz \
  input_R2.fastq.gz \
  result

More on that here:
https://docs.milaboratories.com/mixcr/reference/mixcr-analyze/
After that you will have .clns file for every sample.
Then, to merge all samples you can use mixcr exportClonesOverlap and all .clns files as an input, e.g.:
mixcr exportClonesOverlap \
  --criteria "CDR3|AA|V|J" \
  clonesets/*.clns
  overlapTable.tsv

--criteria allows you to specify the criteria required for clonotypes to be treated as identical. "CDR3|AA|V|J" means that two clonotypes have to share cdr3 amino acid sequence, v gene and j gene.
overlapTable.tsv is exactly what you need.
There are other parameters you can tune (like export only productive clonotypes etc.), more on that here:
https://docs.milaboratories.com/mixcr/reference/mixcr-exportClonesOverlap/#example

Please avoid link only answers, they are useless if the link dies. Summarise the main points of the tutorial in the post. Thanks

Also please disclose if you are in anyway involved with the software you are sharing.

Perfect thank you for the comprehensive update. It is possible Mark Izraelson is employed by MiLaboratories, but we can't be certain.

Hi guys, apologies for initial post. I'm only starting to answer questions on Stackexchange, and will do my best to provide rigorous and useful answers for the community. And yes, I'm involved in MiXCR development for many years so far, so I'm first of all targeting to answer questions related to T and B cell NGS data processing.

@MarkIzraelson no worries, the rules here are a bit stricter than other internet spaces :) And not accusing you of any malfeasance, it's just part of the rules here to disclose affiliation. But it's great to have software developers here answering questions, so we welcome you to keep answering!


text:
Merging two .txt files
I have two .txt files and let's call it 'File A' and 'File B.Here are examples of what file A and file B looks like:
FILE A
#Chr    Start   End Ref Alt Phenotypes
1   1000    1000    A   T   Pheno A
2   2000    2000    T   G   Pheno F, Pheno G
3   3000    3000    C   G   Pheno L

FILE B
#Chr    Start   End Ref Alt Phenotypes
1   1000    1000    A   T   Pheno B
2   2000    2000    T   G   Pheno H, Pheno I
3   3000    3000    C   G   Pheno M
3   3000    3000    C   G   Pheno N

Here's how I want the output to look like:
#Chr    Start   End Ref Alt Phenotypes
1   1000    1000    A   T   Pheno A, Pheno B
2   2000    2000    T   G   Pheno F, Pheno G, Pheno H, Pheno I
3   3000    3000    C   G   Pheno L, Pheno M, Pheno N

Both are tab-separated .txt files and I am trying to merge them together for SNP annotation and would love to hear suggestions on how to proceed this.

This can be also done by Python pandas quite easily, but it you'll need to load each ouput into a dataframe, which is a bit complicated.

Are the #FILE A and #FILEB lines actually part of the file? And is #OUTPUT required in the output? Also, are the files small enough to fit at least one of them in memory? What is the field separator? Is it space(s)? Tabs?

I edited to add actual tabs instead of spaces. Please check and confirm that the files are correct. There's no tab between Pheno F, and Pheno G for example, right?

Whereas wouldn't be as fast as awk, here is a rather dirty R solution:
library(dplyr)
library(tidyr)
library(magrittr)

# File 1

# the "fill" parameter is used since not all the rows have the same number of cols
file_1 <- read.table("file_1.txt", skip = 1, fill = TRUE)

# paste() function with argument collapse is applied to the "phenotype"
# columns of each row (these cols start at position 6 and goes until
# the number of cols: dim(x)[2] gives the number of cols of a table
phenotypes <- apply(file_1[, 6:dim(file_1)[2]], 1, function(x) paste(x, collapse = ""))

# Cols referring to the genomic position, ref and alteration (cols 1:5)
# are combined with the newly created "phenotypes" column
file_1 <- data.frame(file_1[, 1:5], phenotypes)

# colnames are added, could be avoided if dealt with at file reading step
names(file_1) <- c("Chr", "Start", "End", "Ref", "Alt", "Phenotypes")

# File 2

file_2 <- read.table("file_2.txt", skip = 1, fill = TRUE)

phenotypes <-apply(file_2[, 6:dim(file_2)[2]], 1, function(x) paste(x, collapse = ""))

file_2 <- data.frame(file_2[, 1:5], phenotypes)

names(file_2) <- c("Chr", "Start", "End", "Ref", "Alt", "Phenotypes")

# The two tables generated as above are merged using the appropriate
# columns as "anchors" in the "by" argument of full_join()
# The resulting Phenotypes.x and Phenotypes.y cols are concatenated
merged_tables <- full_join(file_1,
                           file_2,
                           by = c("Chr", "Start", "End", "Ref", "Alt"))
merged_tables$Phenotypes <- paste(merged_tables$Phenotypes.x,
                                  merged_tables$Phenotypes.y,
                                  sep = ",")
merged_tables$Phenotypes.x <- NULL
merged_tables$Phenotypes.y <- NULL

# Duplicated rows (in terms of genomic pos, ref and alteration) in any
# of the files result in duplicated rows in the merged table and are
# "handled" by aggregating/summarizing the info at the phenotype cols
# by using paste()
# strsplit() is used to "split" the "phenotypes" on "," in order to
# remove repeating phenotypes resulting from the full_join() call

merged_tables %<>% group_by(Chr, Start, End, Ref, Alt) %>%
      summarize(Phenotypes = paste(Phenotypes, collapse = ",")) %>%
      ungroup() %>%
      mutate(Phenotypes = strsplit(Phenotypes, ","))

merged_tables$Phenotypes <- lapply(merged_tables$Phenotypes, function(x) paste(unique(unlist(x)), collapse = ","))
merged_tables$Phenotypes <- unlist(merged_tables$Phenotypes)

# A tibble: 3 x 6
    Chr Start   End Ref   Alt   Phenotypes
  <int> <int> <int> <fct> <fct> <chr>
1     1  1000  1000 A     T     PhenoA,PhenoB
2     2  2000  2000 T     G     PhenoF,PhenoG,PhenoH,PhenoI
3     3  3000  3000 C     G     PhenoL,PhenoM,PhenoN

# The resulting table (tibble in this case) can be saved as a file
# using write.table() or write.xlsx()

paste() in combination with its collapse parameter is used because of the unequal numbers of columns in the files. fread() of the data.table package could be helpful for reading large files.

This would really benefit from some explanation. I don't have any idea what it's doing (my R-fu is weak). Does it at least produce a final text file as the OP requested? It looks like you're just manipulating the data in R.

I have added quite some comments within the code, explaining each step. The resulting table can easily be exported as a comma or tab separated file via write.table().

Thank you. That's very helpful, have an upvote!

Convert your text files A and B to true, sort-bed-sorted BED and then use bedmap to map their associated IDs:
$ awk -vOFS="\t" '($1!~/#/){ $3+=1; print $1,$2,$3,$4; print $1,$2,$3,$5; print $1,$2,$3,$6; }' A.txt | sort-bed - > A.bed
$ awk -vOFS="\t" '($1!~/#/){ $3+=1; print $1,$2,$3,$4; print $1,$2,$3,$5; print $1,$2,$3,$6; }' B.txt | sort-bed - > B.bed

Then map:
$ bedops -u A.bed B.bed | bedmap --echo --echo-map-id-uniq --delim '\t' --multidelim '\t' - > answer.bed

If you need to put the header back in:
$ bedops -u A.bed B.bed | bedmap --echo --echo-map-id-uniq --delim '\t' --multidelim '\t' - | cat <(echo -e "#OUTPUT\n") <(awk '(NR==2)' A.bed) - > answer.txt

Some work might be needed with awk to put reference and alternate alleles back into their own columns, but this should get you about 90% of the way there.

I just tried this and (assuming I'm using the right executables, I found https://github.com/bedops/bedops) I am getting completely different results to what the OP asked for. The lines are interleaved (some have a residue as the 4th field, others the string "Pheno"), none of the phenotypes are collected, each phenotype has been split into Pheno and N but the two don't sort together so they're all over the place (e.g. one line is 3 3000 3001 C C G Pheno). I don't think this is going to be very helpful to the OP.

Thank you for the suggestion @Alex Reynolds but unfortunately, your commands did not produce the result I wanted but I did learn a thing or two using bedops. Cheers.

As mentioned, you'd need to do some more work with any BED-based approach. An advantage is less memory usage; hash tables in R and other languages tend to use lots of memory. R, in particular, uses a lot of memory. Still, if scalability is not an issue (if your inputs are small) then R is fine. Otherwise, bedops and awk can be useful for a scalable approach.

But this really isn't even close to what the OP was asking for. Could you demonstrate how this can be done using this approach? The output these commands produce is as far from what the OP is asking for as the original input is. What's more, the commands you suggest will mess up the data, mixing words up.

Edit: I could not get bedtools to produce exactly the right output. See note at the end.
First install bedtools. Then, in Bash:
(tail -n +2 FileA.txt; tail -n +2 FileB.txt) | sort -k1,1V -k2,2n -k3,3n - > FileC.txt
head -n 1 FileA.txt > FileD.txt
bedtools merge -i FileC.txt -c 4,5,6 -o distinct,distinct,distinct -delim ", " >> FileD.txt

First line: tail takes lines 2+ of the file, sort sorts the files based on the chromosome, start, and stop coordinates, - is the output from the previous command.
Next line: head takes the first 1 line of FileA (header), saves to FileD.
Last line: bedtools merge merges the intervals that overlap. Option -c tells merge what columns to do operations on when two intervals are merged. Option -o tells merge what to do with those columns; in this case, it will make a comma-delimited list of distinct values. Option -delim defines the delimiter for listing distinct info in columns 4-6.
Output:
#Chr    Start   End Ref Alt Phenotypes
1   999     1001    A   T   Pheno A, Pheno B
2   1999    2001    T   G   Pheno F, Pheno G, Pheno H, Pheno I
3   2999    3001    C   G   Pheno L, Pheno M, Pheno N

PROBLEM: This code will not produce the desired output b/c bedtools expands the interval 1    1000    1000 to 1    999    1001. BED files are used for multi-nucleotide intervals whereas the VCF standard is better suited for single nucleotide intervals and "ref/alt" variant information. Bedtools might handle VCFs with 1-nt intervals appropriately, but but evidently not BED files with 1-nt intervals.

This won't work since the OP's file also has a header and if you sort it, the headers from both files appear as lines in the end. And it also fails since there's no 7th column. I've now edited the question with correctly formatted input files (tabs where they should be). Your answer will fail on those. I think the Phenotype fields won't have tabs in them since they're all the same field. But in any case, this doesn't seem to work here.

Oh man. Yeah, I miscounted the columns, and I wasn't aware sort would take multiple files. AND I posted my answer on a Windows machine, so I (foolishly) didn't test it. I still think bedtools merge is the most simple/elegant way to go. I'll probably implement tail to fix the header problem when I can test it later today.

dammit, bedtools should be the most simple and elegant, indeed. But I've barely used it and both of the answers posted so far get it wrong. Looking forward to what you'll post when you're at a machine where you can test it! If you give me a ping when you do, I'll come and upvote :)

@terdon yeah, it doesn't work. I fixed the problems you mentioned and then found that bedtools "interprets" the intervals. I could write a band-aid for it, but then it becomes as complex as other answers here, so the benefit is lost.

If your files are small enough to fit in memory, here's one way, in Perl:
$ perl -lane 'if(/^#/){$h=$_; next}; $k=join("\t",@F[0..4]); $v=join(" ",@F[5..$#F]); $d{$k} ? ($d{$k} .= ", $v") : ($d{$k} = $v) ; }{ print $h;for(sort keys(%d)){print "$_\t$d{$_}"}  ' file*
#Chr Start End Ref Alt Phenotypes
1 1000 1000 A T Pheno A, Pheno B
2 2000 2000 T G Pheno F, Pheno G, Pheno H, Pheno I
3 3000 3000 C G Pheno L, Pheno M, Pheno N

Since that might be a little unclear, here's the same thing as a commented script and with informative variable names:
#!/usr/bin/env perl

## Iterate over each line of each input file. This is done by the -n switch.
while (<>) {
  ## remove trailing newlines. This is done by -l
  chomp;
  ## If this line starts with a #
  if(/^#/){
    ## Save it as a header. We do this for each header line, so only
    ## the last will be saved. but since the files both have the same header,
    ## that's not a problem.
    $header = $_;
    ## Skip to the next line. This is to avoid processing the headers.
    next
  }
  ## Split the line on whitespace into the @fields array. This is
  ## done by the -a switch which creates the @F array used above.
  @fields = split(/\s+/);
  ## The $key variable holds the concatenated string of the first 5 fields.
  $key=join("\t",@fields[0..4]);
  ## The $values variable holds the rest of the line.
  $values=join(" ",@fields[5..$#fields]);
  ## If we've seen this $key before
  if($data{$key}){
    ## Append the current data to the value stored in the %data hash for this
    ## $key.
    $data{$key} .= ", $values";
  }
  ## If we haven't seen this key before
  else{
    ## Save the values for this key in the %data hash.
    $data{$key}=$values
  }
}
## print the header
print "$header\n";
## Now that all the data have been loaded, iterate over
## the sorted list of keys from the %data hash.
foreach my $key (sort keys(%data)){
  ##And print the output
  print "$key\t$data{$key}\n"
}

If you save that as foo.pl, you can run it like this:
$ perl foo.pl fileA fileB
#Chr Start End Ref Alt Phenotypes
1 1000 1000 A T Pheno A, Pheno B
2 2000 2000 T G Pheno F, Pheno G, Pheno H, Pheno I
3 3000 3000 C G Pheno L, Pheno M, Pheno N

Note that I am assuming the headers #File A and #FileB aren't actually part of your input and that you don't want #OUTPUT included in the output file. If this isn't true, let me know and I can modify this.

Thank you for the reply @terdon. You are correct, #File A and #File B are not part of the header as well as the #OUTPUT. Your perl script works perfectly and thank you for explaining each step. The output is just the way I wanted except that the delimiter is not in tab form (which is my fault for not mentioning it in the post). Is it possible to add this in the script?

@Hadi sure, but what do you mean, exactly? Can you edit your question and show an accurate example of your input file, with tabs? That's why I asked in the comments ;) If you already have tabs in the input, the whole thing is much simpler.

I have edited the original post for accurate examples and yes, my original files use tab as delimiter but output uses single space

@Hadi see updated answer. It should work as you want now.


text:
Visually summarizing phylogenetic networks
After having worked with Bayesian phylogenetic tree inference for some time, I am now trying out a method to infer network phylogenies. To get an idea of what it infers, I'd like to find a way of summarizing the results.
The output of the inference is a NEXUS file containing a stream of networks encoded in eXtended Newick as trees with duplicate node names, like this:
(((((((Kiraman:0.028228,Kui:0.028228)S1:0.0797301)#H4:0.35229)#H2:23.4168)#H1:15.8027,(Kamang:29.3316,#H1:5.45455[&&NHX:gamma=0.36825381342449554])S11:10.3482)S12:1.27331,((Klon:1.16727,((#H4:0.524833[&&NHX:gamma=0.33697036713230577],Kafoa:0.632791)S6:0.0615575)#H3:0.472921[&&NHX:gamma=0.6642547809404831])S2:35.6082,(#H2:22.2107[&&NHX:gamma=0.7544913491697585],((A-Petleng:0.421975,(A-Takalelang:0.168615,(A-Ulaga:0.145472,A-Atimelang:0.145472)S4:0.0231422)S5:0.253361)S7:0.902144,#H3:0.629771)S3:21.3468)S8:14.1045)S9:4.17761)S10:54.5301):0;

￼
There is no formal restriction on the shape of the network, and hybridization edges have a weight 0<γ<1. How do visually or formally summarize a stack of several thousands of these networks, and in particular the backbone tree structure plus a glimpse of the hybridization events they display?
I have played with some tree summarizing approaches and tried to first infer a summary tree and then glue some cross-edges in there, but the way I came up with did not look reasonable.
I am now considering to calculate distances between leaves and using network visualization tools like SplitsTree or NeighborNet to get a visual representation of that, but how shoud I represent the uncertainty of the distances in even a single network (not to mention thousands of them) in the input for these algorithms?
Are there any good approaches to summarize phylogenetic networks of the shape given above? Or to consistently simplify them to a shape that can be easily summarized?

Firstly you treat this as two separate independent calculations.

The traditional approach is simply to produce a consensus phylogeny using your Bayesian method and assess the equivocal bootstraps/Bayesian probability.
Feed your data set (not the trees) into NeighborNet within SplitsTrees. There are statistics that classify the level of departure from a bifurcating phylogeny.

Examine both outputs for orthogonal support for incongruence.
Note Phylogenetic network statistics are a relatively soft area of phylogeny (they are not used as hard evidence).
Networks from Bayesian trees
The final approach (which is what you are asking about), is to feed your Bayesian trees into Dendroscope, this is a GUI that represents the uncertainty between different bifurcating phylogenies as networks and you can use it to understand the incongruence between different trees within your Bayesian MCMC, because the phylogenetic networks directly corresponds to the level of incongruence.
Dendroscope is a nice Java API and is very easy to use.


text:
PDB codes for misfolded proteins?
I'm searching for PDB code examples for "healthy proteins" and misfolded ones, so I can compare the 3D structure.
I'd prefer pairs of a "protein misfolding disease", but could be any. Can I find such a list or how should I search for the particular codes?
I'm particularly interested in ones that also cointain SS-bonds.

Do you mean disulphide bonds between cysteine residues ? These are really rare unless it is chemical modifcation

I'm not sure they are rare, when opening random PDB files, I often come accross disulphide bonds...

Hi @JakeB there is a misunderstanding. Disulphide bonds are incredibly well conserved, even in RNA viruses. Misfolding resulting from a disulphinde bond is rare.

PDB 5E83 sickle cell anaemia.
All details for comparison to wildtype rbc are found in,
Oksenberg, D., Dufu, K., Patel, M.P., Chuang, C., Li, Z., Xu, Q., Silva-Garcia, A., Zhou, C., Hutchaleelaha, A., Patskovska, L., Patskovsky, Y., Almo, S.C., Sinha, U., Metcalf, B.W., Archer, D.R.
GBT440 increases haemoglobin oxygen affinity, reduces sickling and prolongs RBC half-life in a murine model of sickle cell disease.
(2016) Br.J.Haematol. 175: 141-153

Searching "misfolded" on the RCSB PDB turns up 64 hits. Most of these won't be relevant, but there are some interesting entries on tumor associated EGFR and misolded RNA.
In order to experimentally obtain a structure of a protein it generally has to be in a stable conformation, i.e. at low energy. Being soluble also helps. This usually corresponds to the "healthy" form, as misfolding would lead to aggregation.
That said, the structures of certain aggregates have been obtained. For example, see this recent article.
There is also the concept of mutating a protein sequence to fold to a different structure, see for instance here, though I guess you are looking for cases with the exact same sequence.

There are very few real structures of soluble misfolded proteins, for the reasons described by jgreener. An example with a disulfide bond is beta-2-microglobulin, eg PDB entries 2XKS vs 2XKU.
Atomic-resolution structures of amyloid fibrils may be what you're looking for, such as 6GK3 for beta-2-microglobulin amyloid (which has its disulfide intact). A search for "amyloid fibril" at rcsb.org will find many structures. But many amyloid precursor proteins are intrinsically disordered, or are fragments of native proteins.
There's no list that I know of, but here's a recent review: https://www.ncbi.nlm.nih.gov/pubmed/30237470


text:
Total read count of shotgun affected by few highly abundant sequences
I am building statistical models to analyse output from Illumina shotgun sequencing (HiSeq 4000) on stool samples (but RNA-seq data should behave similarily). The raw counts are a statistical sample of the total collection of DNA molecules used as input for the sequencing machine, and I think we can assume that the sequencing platform has limited capacity, such that detected reads are always a small fraction of the the original sequencing library.
Consider two different samples: one control sample and one treatment sample, which is identical apart from a single species (or gene) that has much higher absolute abundance relative what it in the control. As a result, the reads of this species will be much higher in the treatment sample, and also the 'real estate' of the remaining genes in that sample will be decreased. This is simply the nature of relative abundances.
However, size factor estimation that uses the "median ratio method" (http://dx.doi.org/10.1186/gb-2010-11-10-r106) for instance, are based on "experience with real data [...] [that] shows [that] a few highly and differentially expressed genes [or abundant species] may have strong influence on the total read count." (http://dx.doi.org/10.1186/gb-2010-11-10-r106). I am curious to see experimental evidence for this last statement. Or does someone have at least an explanation of how total read count can be affected by the absolute abundance of a species/gene?

That sentence in the Anders and Huber paper is a bit misleadingly worded (I don't think that was their intent). What they mean is that highly expressed and highly differentially expressed genes can out-compete the majority of the other genes, which thereby leads to them appearing to have lower expression if one uses total read counts for normalization (basically, what you said in the second paragraph).
An alternate interpretation of what they wrote (likely the one you understandable understood to have been their meaning) is that some highly expressed genes will inflate the total number of reads a sample actually gets. I have yet to see this in practice, but it certainly COULD happen. There are technical biases in sequencing such that some sequences (often due to length and GC content) are easier to sequence on common machines. A good example of this would be PhiX, which forms clusters much more easily on Illumina machines than real samples. I expect that this is relatively rare in practice, but using a more robust normalization would still help in protecting against it.


text:
Chromosome information for miRNA
Apologies if this question is previously asked. I have a BAM file alignment to hg19 and miRNA gff3 file from mirbase for annotation of BAM file. I want to extract information in the following format (e.g.):
miRNA        chr    chr_start    chr_end    expression_count    sequence
hsa-mir-155  chr2       ...         ....         1020            ATGC....

I know that I can use htseq-count for BAM annotation and for expression counts but how can I get their chromosomal location. The reason for doing this specific way is that I want to show them in a CIRCOS plot.
Just counting is not the goal, their chromosomal location also needs to be extracted. So just subread/htseq-count is not sufficient for extracting genomic location information.
Previously, I had the same idea about processing the sam records. If I want to annotate chromosomal information with miRNA using .gff file, I have to check for each line in .sam records whether their chr_start and chr_end lie in the bandwidth of start and end of any miRNA in .gff file. This process can be very heavy and slow because .sam records are of more than 5 GB. That's why I want to know if there is any extremely fast and generalized way to do it.
I want to avoid processing sam records in python because my experience with extremely heavy data like fastq or sam in python is not good. I would rather prefer awk/bash or any low-level language.
Thanks.

you can do samtools view pathtobamfile, this will convert the file to .sam format and send it to stdout, the columns can then be parsed out which contain the information you are after.

The gff should have the coordinates for each feature.  I'd use featureCounts from subread to count up how many reads hit each feature.


text:
Are there any resources for calculating / representing evolutionary relationships within gene clusters?
I am working with a gene cluster found in several species. The genes involved are pretty much the same, with the odd insertion, but they can be often found in different orientations. So, I can't just do a multiple alignment of the whole cluster. I've done phylogenies of individual genes and I've looked at cluster architecture (indels, rearrangements) and things make sense (eg. limited HGT), but..
...Is there a simpler way to do this? A way to calculate / represent a phylogeny of the whole cluster? Like...some sort of web, maybe? Most of the software that does take rearrangements into consideration is meant for whole genomes, and is unreliable when it comes to fragments about 20-30kb long (like mine).
Thank you for your help!

Can you give us some example data? It's kind of hard to understand what you're looking at. Are you trying to build a phylogeny based on the entire cluster or have you separated each species' cluster into a multifasta file with each gene separate? If you haven't, wouldn't that be a reasonable approach?

If there'd be a way to build a phylogeny based on the entire cluster (that can deal with indels/inversions/rearrangements), I'd love to know it.
I guess I could separate each species' cluster into a multifasta file, but I'm not sure what I could do after that? Could you please explain? It's unfortunate, because I'd lose positional information (eg. let's say species X is closely related to species Y and less so to species Z. In X and Y, the genes are in the order A-B-C-D, in Z they're in the order A-D-C-B), but if it's the only way...

Well, I don't really think the order will make much of a change in the phylogeny. I'd just use the multifasta approach to build the tree and then maybe annotate it with the directionality.

Its easy ... I strongly suspect these are bacteria ...please confirm. The likely application is resistance cassettes again please confirm. Finally this is likely a "genetic island" rather than a single gene ... again please confirm.

@terdon - Thank you. I've never built trees with multiple genes at a time, before. Could you recommend some software for it, please? Until now I've mostly used Seaview and the CIPRES webserver.

@Michael G. - Bacteria, yes. (Silly me, forgot to specify); resistance cassettes, no. And, yes I'm talking of a group of genes rather than a single gene.

How many taxa? 1000s or <100 that is kinda important. I get the idea the fragment is 20-30 kb

@MichaelG. Around 20 taxa.

Oh that's easy. I'll post a full response later/tomorrow if no one else does.

Firstly you wanted a webserver to build a tree, ok, I would use PhyML server. The good thing is that with only 22 taxa the calculations are easily doable.
What I'm assuming is you can map your rearrangements to produce a single homologous gene alignment for all taxa. So you appear to infer you can break your 20-30kb contig into individual aligments of homologous genes which has incorporated the rearrangement and enables easy production of an alignment file. This would leave you with around 20-30 gene alignment files. I hope that is correct, because that is how you get around the rearrangement problem and enables the data to be reassembled for the final tree (below).
Doing the tree is simple, because you can exploit the assumptions used in tree building regarding independence between mutations at each alignment position. The theoretical proof, the data is sound, is complicated. In summary,

build a consensus phylogeny (one tree representative of lots of trees) using two different approaches
combine the data into a single alignment, if the above looks good, or excluding the alignments that are affected recombination.

The theoretical issues with what you want to do are complicated because you need to distinguish bifurcating (clonal) homologues from recombination that could have occured in the reassortment events. The following is a  purists approach, you might not want to do this stuff, but it gives you some idea.
Just before I do that ... I just want to check the alignments assumption is correct? and which bacteria?


text:
Identification of differential genes across 8 groups
I need to identify differential genes across 8 groups. I know this can be done using DESeq2 or EdgeR as these are better suited and equipped.
However, I was thinking if something like this piece of code could be tweaked. Here it works for two groups each with 3 replicates and generates a heatmap. The code is in R
pvals <- apply(dm, 1, function(x) wilcox.test(x~c(rep(0,3), rep(1,9)))$p.value)
heatmap(dm[pvals<0.05,])

Thank you for your inputs and suggestions.

What is your reason not to use DESeq2 or EdgeR? Their models will work much better for RNAseq data than a simple Wilcoxon test.

Welcome to Bioinformatics.SE! At the moment, your question simply asks "what about this?" while ignoring a tremendous amount of scholarship on differential expression analysis. Unless you can more clearly present a specific question or concern, there is a good chance this question will be closed as too broad or unclear.

For a single measurement of some fixed, observable quantity like height or grades, using a Wilcoxon test isn't a bad idea. But RNA-seq is a very complicated system and can't simply be analyzed by something like a Wilcoxon test. Your results will be full of mistakes if you do this. I would recommend that you read some introduction to RNA-seq papers or textbooks and watch some presentations by Wolfgang Huber, Gordon Smyth, and Lior Pachter to see why this is the case


text:
How to handle control samples in CLIP-Seq
I have a CLIP-Seq dataset I'm processing, which includes control samples and no inputs.  This is the second CLIP analysis I've performed to help out users of our genome core facility and the first one only had inputs. I've done a number of ChIP analyses in the past, which I loosely consider to be a similar type of analysis, and they all used inputs, so I've never had to look into the differences in using inputs versus controls when doing peak calling before.  Now that I'm trying to figure out what I need to do, I'm a little lost.
Here's what I've got:

Samples of cross-linked RNA/protein complex that went through a poly-A and then a FLAG pull-down (protein G beads with an antibody).  The complex was eluted, treated with proteinase K, and the RNA was sequenced.
Controls went through the same processing as the samples except the control beads had no antibody

There are replicates in both cases.  (As a side-note, the RNA was treated with RNase III while bound to the beads and the RNA's trimmed off the complexes were also sequenced.  I'm not sure what to do with those samples yet.)
I was going to use one of the CLIP peak callers available on Galaxy, but whichever one I use, all the protocols I've found and the documentation I've looked at use inputs for the peak calling method and not controls.  So what do I need to do differently (if anything) to use the controls as the background in the peak calling?  Or do I need to suggest they also sequence the inputs?
And just so I know (dumb question here)... Is "input" (as the name suggests) the sample as it was before the FLAG pulldown?  Why would one use one or the other (input versus control).

Did you get many reads for your control samples? I would think not, since there was nothing on the beads for them to even stick to. Normally this would be treated as input, i.e., not run through beads at all. What you have is roughly equivalent to using IgG as input in ChIP-seq, which generally works poorly due to few reads. What is the goal of the experiment?

There are roughly 10x fewer reads for the controls. I didn’t design the experiment. I’m just helping them out. I believe they want to know what RNA their protein is binding is binding to and where it’s binding (I.e. what motif).

You'll probably just ignore the control samples then.

I've not handled CLIP-seq before, but we have done iCLIP and in this we had controls, but no input. @DevonRyan is correct, these control samples are similar to IgG pulldowns in ChIP seq.
The way we handled this was to call our peaks on isolated samples. The original peak calling alog from the first CLIP seq paper doesn't require input. Also I believe that clipper doesn't need input (although it will take it). We then removed any peaks from the sample that overlapped a peak in the control.
But most of our analyses we didn't do on peaks - we did lots of things like looking at metagenes, or classes of genes, or particular types of exon etc, and when we did this, because we were averaging over many sites, we could generate the pattern for the control as well as for the input, and often normalise the one to the other.

Is it similar to IgG pulldown though?  The beads had no antibody, thus any reads coming out would be anything that just happened to stick to the raw (protein G) beads.  An IgG pulldown would at least have something.  As I understand it, an IgG antibody is non-specific?  The controls have an order of magnitude fewer reads... So is there a way I can evaluate the effectiveness of the control?

My understanding is that IgG is not non-specific, rather it is specific for the IgG of the host organism from which the primary antibody originated. Thus, in theory it shouldn't bind anything in the experiment, thus the IgG control in a ChIP pulls down two populations of fragments - those bound by a protein that interacts with igG (but shouldn't) and those that interact with the beads. Like CLIP, IgG controls in ChIP have much fewer reads. Despite this, you will almost certainly find that there are peaks in the control - there are definately sequences that interact with the beads.

It seems straightforward that you could call peaks in the control and then exclude those peaks from your analysis of the actual sample - as they would represent likely background signal.


text:
Comparing RNA-Seq and Ribo-Seq
I'm attempting to compare total RNA-seq with Ribo-Seq, to determine if changes in Ribo-Seq are due to changes in transcriptional expression (akin to the analysis performed by Anota2Seq). I am, however, struggling to understand how to accomplish this in DESeq2.
My design is as follows:
cts <- read.delim("Interaction.txt", row.names=1)
coldata <- read.delim("Sample_Sheet.txt", row.names=1)
cts <- round(cts)
dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~Genotype + Type + Genotype:Type)
dds <- DESeq(dds)
resultsNames(dds)

> coldata
                Genotype Type
C_24h_Ribo_1     Control Ribo
C_24h_Ribo_2     Control Ribo
C_24h_Ribo_3     Control Ribo
CHIK_24h_Ribo_1    Viral Ribo
CHIK_24h_Ribo_2    Viral Ribo
CHIK_24h_Ribo_3    Viral Ribo
C_24h_RNA_1      Control  RNA
C_24h_RNA_2      Control  RNA
C_24h_RNA_3      Control  RNA
CHIK_24h_RNA_1     Viral  RNA
CHIK_24h_RNA_2     Viral  RNA
CHIK_24h_RNA_3     Viral  RNA

> resultsNames(dds)
[1] "Intercept"                 "Genotype_Viral_vs_Control" "Type_RNA_vs_Ribo"          "GenotypeViral.TypeRNA"

I do not understand which contrast is the one I need, or how to access the results i'm looking for (or interpret them). Any help?

You would be interested in the interaction term. Did someone else write your code?

Adding to what benn pointed out, at this line:
dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~Genotype + Type + Genotype:Type)

You specified that the counts, can be a explained by Genotype (being Control or Viral) or Type (RNAseq or Ribo) or an interaction effect
What you need is this coefficient, Genotype:Type. To use it, go
res <- results(dds,name="GenotypeViral.TypeRNA")

Based on your setup, you are always comparing RNA to Ribo (see "Type_RNA_vs_Ribo" in resultsNames). If changes to the abundance of a mRNA is similar in RNAseq or Ribo, then the coefficient of the interaction term will be zero and it will be zero in the log2foldchange column in results(res). A positive fold change indicates the log2foldchange between RNA and Ribo is higher in Viral.
If you would like to compare ribo versus and rna-seq, then set RNA as the reference by:
coldata$Type = factor(coldata$Type,levels=c("RNA","Ribo"))
# run deseq again
dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~Genotype + Type + Genotype:Type)
dds <- DESeq(dds)

With this you will compare by default Ribo versus RNA, and the interaction term will tell you how much more enriched or depleted the ribo fraction is, in the Viral samples.
Would recommend checking the vignette out here:
http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#contrasts


text:
Concatenate multiple fasta files into one file per unique entry
I have eight fasta files, one per gene, with 106 fasta entries per file, like this:
gene1.fasta:
>sample1_434564
GATCTCGTAGATGTGA
>sample2_2346764
GATCTCGTAGATGTGA

gene2.fasta:
>sample1_532345
ATGCTGTGAAGTA
>sample2_6040103
ATGCTGTGAAGTA

What I want to do is create one fasta file per sample so that I get all gene sequences for that sample in one file. As you can see, the sample ID is the same before the _ in the fasta headers, but the rest is unique for each fasta file.
Is  there any way to do this with awk or other command-line tools?
Desired output:
Sample1.fasta
>gene1
GATCTCGTAGATGTGA
>gene2
ATGCTGTGAAGTA

Sample2.fasta
>gene1
GATCTCGTAGATGTGA
>gene2
ATGCTGTGAAGTA

What operating system are you using? Do you have access to GNU tools? Is this a Linux machine or something else?

@terdon It's on linux yes

You can do this simply enough in awk:
$ awk -F'_' '{
               split(FILENAME,n,/\./);
               if(/^>/){
                  sample=$1;
                  sub(/>/,"",sample);
                  $0=">"n[1]
               }
               print > sample".fasta"
             }' gene1.fasta gene2.fasta ... geneN.fasta

If I run this on the two files you gave, I get:
$ cat sample1.fasta
>gene1
GATCTCGTAGATGTGA
>gene2
ATGCTGTGAAGTA

$ cat sample2.fasta
>gene1
GATCTCGTAGATGTGA
>gene2
ATGCTGTGAAGTA

Explanation

-F'_' : set the field separator to _ as an easy hack that lets us use $1 for the sample name (sample_NNN, so sample will be $1 and NNN will be $2).
split(FILENAME,n,/\./); : FILENAME always holds the name of the file currently being processed. Here, we are splitting it on . to get the gene's name. This assumes all of your files are named $geneName.fasta. The split creates the array n whose 1st element, n[1] is the gene name.
if(/^>/){ : if this is a fasta header line.
sample=$1 : set the variable sample to the sample's name. See the first note about -F'_'.
sub(/>/,"",sample); : remove the > from the sample variable. We don't want that in the file name.
$0=">"n[1] : set the current line to be > and the gene name.
print > sample".fasta" : print each (modified) line into the appropriate file name.


text:
Can I change the name of file features.tsv to genes.tsv
it is said that the features.tsv from Cell Ranger v3 is analogous to the genes.tsv from Cell Ranger v2. So can I change the file name to genes.tsv for Seurat to read it? I know that Seurat has the update version 3 that can read the features.tsv, and I would like to use the old version. So if I change the file name, did it cause any problem?

You can find an answer here: https://github.com/satijalab/seurat/issues/1196. In brief: Read10X is convenience, you can use the files yourself to generate the matrix but if you only have RNA in your experiment your strategy should work. Did you unzip the file and rename? or rename the zipped file?

Thank you for your reply, and I unzip the file and rename it. Much appreciated!

You can follow the 10X instructions to load your matrix into R:
the following is a direct quote from their webpage:

library(Matrix)
matrix_dir = "/opt/sample345/outs/filtered_feature_bc_matrix/"
barcode.path <- paste0(matrix_dir, "barcodes.tsv.gz")
features.path <- paste0(matrix_dir, "features.tsv.gz")
matrix.path <- paste0(matrix_dir, "matrix.mtx.gz")
mat <- readMM(file = matrix.path)
feature.names = read.delim(features.path,
                           header = FALSE,
                           stringsAsFactors = FALSE)
barcode.names = read.delim(barcode.path,
                           header = FALSE,
                           stringsAsFactors = FALSE)
colnames(mat) = barcode.names$V1
rownames(mat) = feature.names$V1

You can then create a seurat object using mat
object <- CreateSeuratObject(mat, project="")


text:
should I continue the iteration of psiblast if the blast step produced only below threshold sequences?
I used a PSIBLAST (position specific scoring) to detect the similarities between specific proteins and the database of Clostridium.
The results did not reveal any similarities above the e value threshold of 0.001.
Could I continue iterations of the second and third round of PSIBLAST?


text:
Mapping probe ID to gene name
I have set of modules from various condition  ,with first column as Probe_ID name and I have an annotation file with Probe_ID and Gene_name.
I have to map the Probe_ID to gene name in each file and write it back.
So what im doing is as such
filelist = list.files(pattern = ".*.txt")

I get list of files
Then those list of files into list of dataframe
datalist = lapply(filelist, function(x)read.table(x, header=T))

Now i read into the annotation file
ANNOTATION_FILE <- read_delim("/run/media/punit/data1/Dev_bio_info/ANNOTATION_FILE.txt",
    "\t", escape_double = FALSE, trim_ws = TRUE)

The normal way i do is
I perform inner_join(ANNOTATION_FILE,my_file) for single file but as i have like 22 files in the module .I want to do it in loop and im not sure how to proceed further.
My data frame example
This is how my module files looks like .
dput(head(module_ivory[c(1:10)]))
structure(list(Gene = c("1556672_a_at", "1557803_at", "1558080_s_at",
"1566465_at", "213089_at", "213229_at"), GSM1304905 = c(5.977158746759,
6.83852532860115, 5.93931007942615, 6.11926533584282, 9.92182534751875,
10.2576150019626), GSM1304906 = c(7.51622774902859, 7.71268478229187,
7.32360342148457, 6.8928342169826, 9.12807784837542, 10.6767090567628
), GSM1304907 = c(6.58069278884036, 7.86890363253932, 6.9901005950279,
6.08614452085924, 10.1220627669628, 10.9203125158285), GSM1304908 = c(6.60119000844375,
6.71735654004629, 7.43132441110608, 5.48084261893207, 9.81689665710992,
10.7649566455486), GSM1304909 = c(5.95049879217177, 6.132975671606,
7.0902951305213, 5.59408027736444, 9.08474568290096, 9.84784203563294
), GSM1304910 = c(5.83616005353191, 5.94947517647125, 6.61469427098759,
4.31285839654272, 9.44446310110043, 9.99021299504538), GSM1304911 = c(5.92373056182016,
6.14502105445603, 6.2333758614231, 5.46455403294761, 9.55447963883665,
10.1810888944642), GSM1304912 = c(6.43107367158601, 6.36340779490876,
7.18374247647159, 5.71508717609646, 9.61355221684561, 10.3266656836978
), GSM1304913 = c(6.4010171640209, 7.31600259413326, 6.16969520557938,
5.40234538288784, 8.89481478840621, 10.290064262843)), class = c("tbl_df",
"tbl", "data.frame"), row.names = c(NA, -6L))

My annotation file
dput(head(ANNOTATION_FILE))
structure(list(Gene = c("1007_s_at", "1053_at", "117_at", "121_at",
"1255_g_at", "1294_at"), Gene_Symbol = c("DDR1 /// MIR4640",
"RFC2", "HSPA6", "PAX8", "GUCA1A", "MIR5193 /// UBA7")), class = c("tbl_df",
"tbl", "data.frame"), row.names = c(NA, -6L))

I would like to do inner join on each file and write them back into the file .
Any help or suggestion would be really appreciated

You mention that you tried a loop , could you show us the code you used and why it didn't work? Also I'm surprised you needed inner_join. Haven't the annotation for your microarray been helpful to map your IDs. Also probably it is worth to first add the annotation you want with before splitting into files (you seem to be using WGCNA or some sort of clustering algorithm, so you could add the gene name before splitting them into 22 files).

To end up with a single table made by merging all your tables, you can give a try with reduce() of purrr:  The first argument .x would be the list of tables to merge (in your case the annotation file plus all the other files) and the .y argument would be your function, i.e. inner_join(). reduce() would merge all the tables within the given list, left to right or right to left. ... would enable you to specify column(s) to use as "anchors" while merging.
If you would prefer to keep your tables separately then lapply() is the answer: lapply(list_of_files, function(x) inner_join(annotation_file, x)). This would give a list of merged tables.

just a doubt is it going to merge all the tables "would merge all the tables within the given list," as you wrote this .Because i want to keep all the individual files separately

reduce() would merge everything and would provide a single table. If you prefer to have separate (merged) files, lapply() would help. I have edited my answer to reflect these.

"If you prefer to have separate (merged) files" yes this is what i would like to .I will give it a try

I ran to this error "Error: by required, because the data sources have no common variables
Call rlang::last_error() to see a backtrace " even though I have a common column name called "Gene"

Works for me with the limited data you provided in your post (except for getting 0 x 11 tables as there are no shared gene identifiers in your annotations table and your data table, hence inner_join() returning 0 rows). Please edit your original post as to include the code you used.

okay but if i do this A <- inner_join(ANNOTATION_FILE,module_darkgreen)
Joining, by = "Gene" i get a file mapped .Not sure what is the issue

im using the same code as given here


text:
How to find common sequences among 6 multi-fasta files
I have 6 multi-fasta files, every of them contains ca 1500 sequences like that:
>Haladaptatus sp. R4
MQATVRDLNGEDADTVDLPDVFETTVRTDLIKRAVLAAQANRKQDYGTDPHAGMRTSAES
PGSGRGMAHVPQTNGRGARVPFTVGGRVAHPPKAEKDRSRSINKKERKLAVRSAIAATTD
AERVSERGHRFDEDTELPLVVSDDFEDLVKTQEVVSFLEAVGIDADIARAEDNKKVRAGR
GTTRGRKYKTPKSILFVTSEEPSRAARNLAGADVATAREVNTEDLAPGTQAGRLTVWTES
ALEEVADR
>Methanosphaera sp. WGK6
MAKVNVYSLKGDITEEIELPEIFEEEYRPDVIKRAVISTQTARIQPWGANPMAGKRTTAE
SFGSGRGAAMVPRVKGSSKAAFVPQAIGGRKAHPPRVNTIYHEKINKKERILAIRSAIAA
TANKEIVEQRGHAVANLEQVPFVVDDELETIKTTKETREIFKDLGIMDDILRAKKGRKIK
SGKGKLRGRKYRTPKGPLVVVGNDRGISLGARNHAGVEVVEVNNINAELLAPGTHAGRLT
IYTKSAVEKLADLFQQNRS

and so on.
I need to create a new result file containing sequences that are common in these multi-fasta files.
Example:
If sequence
>Haladaptatus sp. R4
    MQATVRDLNGEDADTVDLPDVFETTVRTDLIKRAVLAAQANRKQDYGTDPHAGMRTSAES
    PGSGRGMAHVPQTNGRGARVPFTVGGRVAHPPKAEKDRSRSINKKERKLAVRSAIAATTD
    AERVSERGHRFDEDTELPLVVSDDFEDLVKTQEVVSFLEAVGIDADIARAEDNKKVRAGR
    GTTRGRKYKTPKSILFVTSEEPSRAARNLAGADVATAREVNTEDLAPGTQAGRLTVWTES
    ALEEVADR

occures in 1 and 2 and 3 and 4 and 5 and 6 file it shoud be written to the result file.
May anoyone help?
Edited: this is my Python script which should work on 2 files (not 6 for now to make it simplier):
for z in f1:
if z[0] == '>':
    znaleziony = False
    gatunek1 = z[1:len(z)]
    #print(gatunek1)
    for x in g1:
        #print(x)
        if x[0] == '>':

            gatunek2 = x[1:len(x)]
            if gatunek1 == gatunek2:
                wynik.write('>')
                wynik.write(gatunek1)
                znaleziony = True
                continue

elif znaleziony == True:
    wynik.write(z)

wynik.close()
Unfortunately it repeats some identifiers of sequence and does compare only identifiers (not whole sequence).

Are you looking for duplicates?

Welcome to the site. What approach have you tried to do this? In which step are you unsure how to continue? Could you post your code so that we can advice how to continue?

@MichaelG. I woudnt say so because there shouldn be any duplicates among one separate file

@Ilrs I tried to run a script in Python. I`ve begun with finding common sequences between two files (not six for now) to make it simplier. It finds common sequences but also writes only names of several sequences many times in the file like this: >Natrialbaceae

Natrialbaceae
Natrialbaceae
MEATVRDLDGSDAGSTELPAVFETTYRPDLIARAVRVAQANRKQDYGADEFAGMRTPAES
FGSGRGMAHVPRQEGRGRRVPQTIKGRKAHPPKAEKDQSESINTKEKKLAVRSAIAATTD
AELVADRGHQFDDDAEIPVVVSDEFEDLVKTKEVVEFLEAAGLEADVERADEGRSIRSGR
GKTRGRKYKQPKSILFVTSSESGPSRAARNLAGADVTTAAEVNAEDLAPGAQAGRLTVWT
ESALEEVADR

`for z in f1:
if z[0] == '>':
znaleziony = False
gatunek1 = z[1:len(z)]
#print(gatunek1)
for x in g1:
#print(x)
if x[0] == '>':
gatunek2 = x[1:len(x)]
if gatunek1 == gatunek2:
 wynik.write('>')
 wynik.write(gatunek1)
 znaleziony = True
 continue

elif znaleziony == True:
wynik.write(z)
wynik.close()`

Please edit your post and add the code there - comments don't allow for good code formatting and that makes your code difficult to understand.

Based on your question it is not entirely clear if you are looking for common header + sequence or also sequences that are identical too each other but have been mislabeled. "Unfortunately it repeats some identifiers of sequence and does compare only identifiers (not whole sequence)." I would assume that identical identifiers have the same sequence. Can you elaborate on what exactly you are looking to extract?

Hello @Patrycja, your code is not well indented. Please edit the post and indent the code properly as Python blocks are identified by indentation and bad indentation makes the code impossible to read.

Hi @Patryvia if I get.a.minute I'll do it in python pandas. The duplicated function will make short work of the problem

If both the ID lines and the sequences are duplicated, you can do this using the FastaToTbl and TblToFasta scripts I have posted before and just normal Linux utilities:
FastaToTbl file1.fa file2.fa ... file6.fa |
    sort | uniq -dc | awk '$1==6' | sed 's/^  *6 //' | TblToFasta > dupes.fa

To understand this, let's look at the output of each step. The FastaToTbl script
just converts fasta to a tab-delimited format where the ID is the first field and the sequence is the second, all on one line:
$ FastaToTbl file*fa | head -n1
Haladaptatus sp. R4 MQATVRDLNGEDADTVDLPDVFETTVRTDLIKRAVLAAQANRKQDYGTDPHAGMRTSAESPGSGRGMAHVPQTNGRGARVPFTVGGRVAHPPKAEKDRSRSINKKERKLAVRSAIAATTDAERVSERGHRFDEDTELPLVVSDDFEDLVKTQEVVSFLEAVGIDADIARAEDNKKVRAGRGTTRGRKYKTPKSILFVTSEEPSRAARNLAGADVATAREVNTEDLAPGTQAGRLTVWTESALEEVADR

Then, we pass this through sort and uniq -dc which prints out the number of times a line was repeated (only for duplicated lines):
$ FastaToTbl file*fa | sort | uniq -dc
      4 some other sequence ALASOOPAKASKAPPPQQCDDREETAYSSLDAS
      6 Haladaptatus sp. R4 MQATVRDLNGEDADTVDLPDVFETTVRTDLIKRAVLAAQANRKQDYGTDPHAGMRTSAESPGSGRGMAHVPQTNGRGARVPFTVGGRVAHPPKAEKDRSRSINKKERKLAVRSAIAATTDAERVSERGHRFDEDTELPLVVSDDFEDLVKTQEVVSFLEAVGIDADIARAEDNKKVRAGRGTTRGRKYKTPKSILFVTSEEPSRAARNLAGADVATAREVNTEDLAPGTQAGRLTVWTESALEEVADR

Since we only want those that were present in all 6 files, we can use awk telling it to only print lines whose first field was 6:
$ FastaToTbl file*fa | sort | uniq -dc | awk '$1==6'
      6 Haladaptatus sp. R4 MQATVRDLNGEDADTVDLPDVFETTVRTDLIKRAVLAAQANRKQDYGTDPHAGMRTSAESPGSGRGMAHVPQTNGRGARVPFTVGGRVAHPPKAEKDRSRSINKKERKLAVRSAIAATTDAERVSERGHRFDEDTELPLVVSDDFEDLVKTQEVVSFLEAVGIDADIARAEDNKKVRAGRGTTRGRKYKTPKSILFVTSEEPSRAARNLAGADVATAREVNTEDLAPGTQAGRLTVWTESALEEVADR

But then we don't want the 6 as part of the final output, so we need to remove it:
$ FastaToTbl file*fa | sort | uniq -dc | awk '$1==6' | sed 's/^  *6 //'
Haladaptatus sp. R4 MQATVRDLNGEDADTVDLPDVFETTVRTDLIKRAVLAAQANRKQDYGTDPHAGMRTSAESPGSGRGMAHVPQTNGRGARVPFTVGGRVAHPPKAEKDRSRSINKKERKLAVRSAIAATTDAERVSERGHRFDEDTELPLVVSDDFEDLVKTQEVVSFLEAVGIDADIARAEDNKKVRAGRGTTRGRKYKTPKSILFVTSEEPSRAARNLAGADVATAREVNTEDLAPGTQAGRLTVWTESALEEVADR

Finally, convert to fasta again:
$ FastaToTbl file*fa | sort | uniq -dc | awk '$1==6' | sed 's/^  *6 //' | TblToFasta
>Haladaptatus sp. R4
MQATVRDLNGEDADTVDLPDVFETTVRTDLIKRAVLAAQANRKQDYGTDPHAGMRTSAES
PGSGRGMAHVPQTNGRGARVPFTVGGRVAHPPKAEKDRSRSINKKERKLAVRSAIAATTD
AERVSERGHRFDEDTELPLVVSDDFEDLVKTQEVVSFLEAVGIDADIARAEDNKKVRAGR
GTTRGRKYKTPKSILFVTSEEPSRAARNLAGADVATAREVNTEDLAPGTQAGRLTVWTES
ALEEVADR

Alternatively, if your files are small enough to be loaded into memory, you could do this instead:
FastaToTbl file*fa |
    awk '{a[$0]++}END{for(s in a){if(a[s]==6){print s}}}' |
        TblToFasta > dupes.fa

Or you could just do the whole thing in a single awk operation:
awk '{
        if(/^>/){
            a[seq]++;
            seq=$0;
        }
        else{
            seq=seq"\n"$0
        }
      }
      END{
        a[seq]++;
        for(s in a){
            if(a[s]==6){
                print s
            }
        }
       }' file*.fa

For pandas you simply load the fasta names and fasta sequences into separate arrays. The "fasta names" array becomes the index in a pandas Series, whilst the sequence are the values. The dataframe will correctly re-associated the names with each sequence. You load all the files into a single dataframe and use the duplicated() function then write out.
The easiest way to do this is concatenate all your files into one via shell ...
cat file*.fasta > concat_file.fa

Python script (biopython + pandas) ...
from Bio import SeqIO
import pandas as pd

path = "/Users/username/location/concat_file.fa"
outpath = "/Users/username/Desktop/duplicates.fasta"
fasta_name = []
fasta_seq = []

for record in SeqIO.parse(path, "fasta"):
    fasta_name.append('>' + record.id)
    fasta_seq.append(record.seq)
df = pd.Series(fasta_seq, index = fasta_name)
df_dups = df.index.duplicated()
mydups = df[df_dups]
mydups.to_csv(outpath, sep="\n", index=True, header=False)

...
I checked it on a ~4000 sequence fasta file of unique sequence ids (dengue), duplicated the last sequence id + sequence and the only sequence output was the very last sequence.


text:
how to merge more than two sample in Seurat?
I would like to merge more than two sample in the Seurat, and the mergeseurat can only merge two sample. So what should I do now. The screenshot is my script.

You have fed arguments to the MergeSeurat() function that it does not expect. In terms of objects, MergeSeurat() accepts only 2 arguments object1 and object2, please do ?MergeSeurat and see yourself.
Seurat is probably the best documented single cell package, your exact question has its own dedicated vignette: https://satijalab.org/seurat/v3.0/merge_vignette.html
As the vignette above is for Seurat v3 and in case you would like to use Seurat v2:
You can create a named list of Seurat objects to be merged and then use reduce() of the purrr package.
# create a named list of seurat objects to be merged
# code not shown

# optional but probably a good idea
# rename cells using object names as prefix
for (i in names(seurat_object_list)) {
  seurat_object_list[[i]] <- RenameCells(seurat_object_list[[i]],
                                         add.cell.id = i)
}

# merge all the objects in the list
merged_combined <- reduce(seurat_object_list,
                          MergeSeurat,
                          do.normalize = FALSE)

Hi haci, I have tried, and it got the error "Error in .f(out, .x[[i]], ...) : trying to get slot "raw.data" from an object of a basic class ("character") with no slots". Do you have any idea about this error? Much appreciated!

You are feeding a character vector to your function (probably MergeSeurat() used within reduce() ) and it complains because it is looking for the "raw.data" slot of your object. So you should feed a proper Seurat object to your function. You can use class() to check the nature of your object.

I have check it. And it shows " [1] "seurat" ". I don't know whether it is correct.

What are the outputs of these commands? 1) lapply(seurat_object_list, class), 2) lapply(seurat_object_list, function(x) dim(x@raw.data))

I got "[[1]][1] "character"[[2]][1] "character"[[3]][1] "character"[[4]][1] "character" on your first command. And got the error "Error in FUN(X[[i]], ...) : trying to get slot "raw.data" from an object of a basic class ("character") with no slots" in the second command. Much appreciated.

These prove that your object list does not contain seurat objects but just character vectors which obviously does not have raw.data slots. Please double check the way you create your object list. Should be sth like object_list <- list(name1 = object1, name2 = object2)

The MergeSeurat command is from Seurat v2.4 and only accepts two objects as parameters. If you need to merge more than one you can first merge two, then merge the combined object with the third and so on.
In the newer Seurat v3.0 this is replaced by the merge command that can have a named list of Seurat objects as input
# merge two objects
merge(x = pbmc_small, y = pbmc_small)
# to merge more than two objects, pass one to x and a list of objects to y
merge(x = pbmc_small, y = c(pbmc_small, pbmc_small))

If you want to perform batch correction at the same time as merging the data you should have a look at the IntegrateData functionality of Seurat v3.0, explained in the corresponding vignettes:
https://satijalab.org/seurat/v3.0/immune_alignment.html
https://satijalab.org/seurat/v3.0/integration.html
A similar workflow for Seurat v2.4 can be cound here:
https://satijalab.org/seurat/v2.4/immune_alignment.html

As of 2023 (or Seurat >4 or >3?), simply use merge and write
combined.obj <- merge(x = ls.Seurat[[1]], y = ls.Seurat[2:l(ls.Seurat)], add.cell.ids = vec_sample_IDs)

Since y can take a list.

For now, Seurat does not have an option to merge more than two samples at a time. You can write a function to run the MergeSeurat function in a loop. Following issue in Github has an example function for multiple samples.
https://github.com/satijalab/seurat/issues/432


text:
Convert paired-end BAM into a single-end BAM and keep all the reads
I want to call peaks using MACS2 on my paired-end ATAC-seq BAM and I am interested in cutting sites. So I make use of the -BAM option (opposed to the -BAMPE) option. However this throws away all the mates, and all the reads of which the mate did not align (that's more than 50%!!!). My idea was to align paired-end, and then remove all the information from the BAM that it is paired-end. So my question is how do I:
Convert paired-end BAM into a single-end BAM and keep all the reads
I made a start based on the SAM-format specification , but when I run the result of this script through markduplicates I get the error:
ERROR: Record 1, Read name SRR6245286.60336908., MRNM should not be set for unpaired read.

I can not find what this MRNM is...
This is my code I used to generate the faulty 'single-end' bam, but I am open to suggestions (e.g. bash, preferably no R).
import pysam

paired =         1
proper_pair =    2
mate_unmapped =  8
mate_reverse =   32
first_in_pair =  64
second_in_pair = 128

bam_in = pysam.AlignmentFile("/vol/atacpipeline/bwa/GSM2837489-danRer7.samtools-coordinate.bam", "rb")
bam_out = pysam.AlignmentFile("/vol/atacpipeline/bwa/GSM2837489-danRer7_SE.samtools-coordinate.bam", "wb", template=bam_in)
for line in bam_in:
    if line.flag & second_in_pair:
        line.pos += line.template_length
    line.next_reference_id = 0
    line.next_reference_start = 0
    line.flag &= ~(paired + proper_pair + mate_unmapped + mate_reverse + first_in_pair + second_in_pair)

    bam_out.write(line)

That error is probably from Picard, making it irrelevant. Your script should work fine for macs2.

@DevonRyan Yep, my bad. Everything works just fine as is with macs2. Still weird that picard does not accept the file, but that is another problem.

Picard is incredibly picky :(

Like Devon pointed out, most likely you should sort out whether the files have been marked for duplicates correctly.
You can also use samtools rmdup too
samtools rmdup <bamfile> <output>

Also, I find it very odd that 50% of the reads are lost when you keep only properly paired reads.
I would suggest removing the duplicates and then converting the bam file to a bed file. Macs2 can take in the bed format. Any extremely crude python script that converts a bam into a bed:
import pysam
samfile = pysam.AlignmentFile("test.bam", "rb")
open("test.bed", 'w').close()
f = open("test.bed", "a")
STRAND = ["+","-"]
for read in samfile.fetch():
    STR = STRAND[int(read.is_reverse)]
    BED =[read.reference_name,read.pos,read.reference_end,".",read.mapq,STR]
    f.write('\t'.join([str(i) for i in BED]))
    f.write('\n')
f.close()

You can also try the conversion with bedtools or definitely refine that script above with pandas etc.

The problem with macs2 if you look for cutting sites is that they throw away all the mates (which is 50% of you reads) + only keeps properly aligned. So that is at least 50%. I'll take a look at changing it into bed and see how that works. Thanks

MRNM stands for "Mate reference index". So Picard found something in the RNEXT field which should be set only for paired-end reads but the rest of the file looks like single-end.
The problematic line in your code is:
line.next_reference_id = 0

This sets the RNEXT SAM field to whatever Pysam stores as a reference with index 0 (next_reference_id). This is not the same as reference name (next_reference_name). In single-end SAM file this field should be set to *. You cannot set * with next_reference_id but you can set it with next_reference_name:
line.next_reference_name = "*"

I would also recommend to change next_reference_start line to:
line.next_reference_start = -1

just to be sure Pysam handles this field as not set (I tried to put there 0 Pysam changed it to 1 for some reason). I would also add line:
line.template_length = -1

because this field should also be empty in single-end SAM.
This worked in my case. If Picard keeps complaining you can add VALIDATION_STRINGENCY=LENIENT or SILENT to your Picard command to suppress exiting of the command at errors.


text:
Fix dash/empty ALT alleles for deletions in VCF
I have a VCF file that has deletions specified like this, with the REF allele containing only the deleted portion and the ALT allele just saying -:
chr1 101 . AAA -

I would like to combine this with another VCF in which deletions are specified like this, with REF and ALT both extended by 1bp to the left to ensure neither is empty:
chr1 100 . GAAA G

Is there a tool to convert from one of these representations to the other? I have already tried vt normalize and bcftools norm but neither of them seem to do this.
I don't really mind which representation I end up with although I think the second is the one more consistent with the spec.

Cross-posted on biostars: https://www.biostars.org/p/394866/

OP, the first of the two representations is so outdated that most tools might error out. You're better off insisting on converting to format #2

@RamRS I agree. Would definitely prefer a way to convert format 1 to format 2 than vice versa. The other way round could be done with a simple awk solution though, since it doesn't require knowledge of the reference sequence.

According to the specs the ALT is not allowed to contain -. While technically supported, 'N' as the alternative allele is sometimes used unnecessarily and this script can also fix those cases if one exchanges the two "-" within the fix_record function for "N".
Here is a little python (>=3.6) script that is able to fix it. It makes use of the pysam module. You need a reference file and this fasta file must be indexed by samtools faidx reference.fa. Please note that it doesn't support multiple values in the ALT column for records that contain a -.
Usage:
python fixvcf.py -f reference.fasta -o out.vcf in.vcf

The script:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import argparse
import sys

import pysam

class Application:
    def __init__(self):
        self.args = Application.get_args()
        self.vcf = pysam.VariantFile(sys.stdin)
        self.fasta = pysam.FastaFile(self.args.fasta)

    @staticmethod
    def get_args() -> argparse.Namespace:
        """
        Read command line arguments and set input and output sources.

        :return: parsed command line arguments
        """
        parser = argparse.ArgumentParser(prog="fixvcf.py", description="Fixes deletion represented by -")
        parser.add_argument('--version', action='version', version='%(prog)s 0.1')

        parser.add_argument("input", help="input file, use - to read from stdin")
        parser.add_argument("-f", "--fasta", help="fasta reference file", required=True)
        parser.add_argument("-o", "--output", help="output file")

        if len(sys.argv) == 1:
            parser.print_help(sys.stderr)
            sys.exit(1)

        args = parser.parse_args()

        if args.output:
            sys.stdout = open(args.output, "w")

        if args.input != "-":
            sys.stdin = open(args.input, "r")

        return args

    def fix_record(self, record: pysam.libcbcf.VariantRecord) -> pysam.libcbcf.VariantRecord:
        """
        If the ALT column of the records contains a `-`, the variant position is shift one position to the left.
        For reaching this the reference base for the position is determined and prepended to the current REF bases.
        The ALT value `-` is replaced with the determined ref base on position before. In case there are multiple ALT
        values the determined ref base is prepended to them.

        :param record: vcf file record
        :return: fixed vcf record
        """
        if "-" in record.alts:
            ref = self.fasta.fetch(region=f"{record.chrom}:{record.pos - 1}-{record.pos - 1}")
            record.pos = record.pos - 1
            record.ref = f"{ref}{record.ref}"
            record.alts = tuple(map(lambda alt: ref if alt == "-" else f"{ref}{alt}", record.alts))

        return record

    def start(self):
        """
        Print the header of the vcf file. Afterwards the fix_record method is applied to each record in the vcf file.

        For more information about python's map() see: https://docs.python.org/3.6/library/functions.html#map

        For more information about the unpacking mechanism using `*`see:
        https://docs.python.org/3/tutorial/controlflow.html#tut-unpacking-arguments
        """
        print(self.vcf.header, end="")
        print(*map(self.fix_record, self.vcf), sep="", end="")

if __name__ == "__main__":
    app = Application()
    app.start()

A lot of the code above is boilerplate for running a script on the command line. The critical code for the transformation is in the fix_record function. The relevant code is below, expanded, with line numbers and a description of each command.
[1]   ref = self.fasta.fetch(
          region=f"{record.chrom}:{record.pos - 1}-{record.pos - 1}"
      )
[2]   record.pos = record.pos - 1
[3]   record.ref = f"{ref}{record.ref}"
[4]   record.alts = tuple(
          map(
              lambda alt: ref if alt == "-" else f"{ref}{alt}",
              record.alts
          )
      )

This grabs the base preceding the inserted/deleted sequence.
This corrects the position of the variant to account for prepending another nucleotide.
This defines the new corrected REF allele: concatenate the base grabbed in [1] with the original REF allele
Apply an anonymous (lambda) function to each ALT allele; the function will replace - with the base grabbed in [1], or preprend the base to the allele if it is not -

+1, thanks for sharing a working script. It might be even more helpful for innocent bystanders to emphasize and explain the code that does the actual transformation.

Done by adding type hints and doc strings :)

Sorry, I should have been more clear. My edit is more of what I had in mind. :-)

Thanks finswimmer for the script. It's been very helpful to me. I would like to extend it for the case of REF alleles with "-" as in: chr1 6029086 - A to: chr1 6029086 CA C I have made the attempt to do it by extend this solution for other cases but couldn't find the right way to do it for all of them. For example in: chr1 1684347 - CCT I would change it to (REF would change to "C"(the position grabbed)+ALT and ALT would change to the position grabbed("C") which worked in the previous example): chr1 1684347 CCCT C But it's wrong cause actually should be like: chr1 1684347 C CCCT How

Generator version of the python script. It won't consume your memmory at all. #print(*map(self.fix_record, self.vcf), sep="", end="") for line in self.vcf: print(self.fix_record(line), sep="", end="")

Probably it is better to use dedicated tools, however, here is a quick R workaround. It is relatively straightforward to remove the extra base (first letter of column 4) in the second type of output:
# just a workaround to encolose the second file in a data frame
# you can just read your file as a table using fread(), read.csv(), ...
test <- data.frame(
  col1 = "chr1",
  col2 = 100,
  col3 = ".",
  col4 = "GAAA",
  col5 = "G",
  stringsAsFactors = FALSE
)

> test
  col1 col2 col3 col4 col5
1 chr1  100    . GAAA    G

# remove the first letter in column 4
test$col4 <- substring(test$col4, 2)
# converting column 4 to "-"
test$col5 <- "-"

> test
  col1 col2 col3 col4 col5
1 chr1  100    .  AAA    -

EDIT:
Here is how it can be done with awk for speed and memory concerns: The logic is the same as above, remove the first letter in column 4 and then change column 5 to "-".
awk '{$4 = substr($4,2); print}' infile > outfile
awk '{$5 = "-"; print}' outfile > outfile2

infile
chr1 100 . GAAA G

outfile
chr1 100 . AAA G

outfile2
chr1 100 . AAA -

Thanks but my real VCF is several GB in size so I don't think this will work.

I have edited my answer so as to include a solution making use of awk.

One more thing. If empty ALT allele is not "-", you can change it such as "*", as follows.
    if "*" in record.alts:
        ref = self.fasta.fetch(reference=record.chrom, start=record.pos - 1, end=record.pos)
        record.pos = record.pos - 1
        record.ref = f"{ref}{record.ref}"
        record.alts = tuple(map(lambda alt: ref if alt == "*" else f"{ref}{alt}", record.alts))

Hi Namshin, thanks for contributing, but this is not the way we usually answer questions on this website. You should not assume that a reader reads the whole thread. If it's an extension for a different answer, meaning which one and also explain what is the benefit of your extension. Note you can [[edit]] your answer.

This is my first post and I will add some explanations. When I tried to use this script, print(*map(self.fix_record, self.vcf), sep="", end=""), it will increase memory usage because it will load everything into memory and then start to apply fix_record function. So, I just used two lines for line in self.vcf: print(self.fix_record(line), sep="", end=""). Because pysam uses generator by default, if we process one by one by pysam's VariantFile function, it will consume very little memory. In my case, the VCF file was more than 200GB, so I had to fix this script for my VCF file.
Second my VCF file used "*" for empty ALT, so I fixed two lines (1) if "*" in record.alts: (2) record.alts = tuple(map(lambda alt: ref if alt == "*" else f"{ref}{alt}", record.alts))
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import argparse
import sys

import pysam

class Application:
    def __init__(self):
        self.args = Application.get_args()
        self.vcf = pysam.VariantFile(sys.stdin)
        self.fasta = pysam.FastaFile(self.args.fasta)

    @staticmethod
    def get_args() -> argparse.Namespace:
        """
        Read command line arguments and set input and output sources.

        :return: parsed command line arguments
        """

        parser = argparse.ArgumentParser(prog="fixvcf.py", description="Fixes deletion represented by -")
        parser.add_argument('--version', action='version', version='%(prog)s 0.1')

        parser.add_argument("input", help="input file, use - to read from stdin")
        parser.add_argument("-f", "--fasta", help="fasta reference file", required=True)
        parser.add_argument("-o", "--output", help="output file")

        if len(sys.argv) == 1:
            parser.print_help(sys.stderr)
            sys.exit(1)

        args = parser.parse_args()

        if args.output:
            sys.stdout = open(args.output, "w")

        if args.input != "-":
            sys.stdin = open(args.input, "r")

        return args

    def fix_record(self, record: pysam.libcbcf.VariantRecord) ->     pysam.libcbcf.VariantRecord:
        """
        If the ALT column of the records contains a `-`, the variant position is shift one position to the left.
        For reaching this the reference base for the position is determined and prepended to the current REF bases.
        The ALT value `-` is replaced with the determined ref base on position before. In case there are multiple ALT
        values the determined ref base is prepended to them.

        :param record: vcf file record
        :return: fixed vcf record
        """

        if "*" in record.alts:
            ref = self.fasta.fetch(reference=record.chrom, start=record.pos - 1, end=record.pos)

            record.pos = record.pos - 1
            record.ref = f"{ref}{record.ref}"
            record.alts = tuple(map(lambda alt: ref if alt == "*" else f"{ref}{alt}", record.alts))

        return record

    def start(self):
        """
        Print the header of the vcf file. Afterwards the fix_record method is applied to each record in the vcf file.

        For more information about python's map() see:
https://docs.python.org/3.6/library/functions.html#map

        For more information about the unpacking mechanism using `*`see:
        https://docs.python.org/3/tutorial/controlflow.html#tut-unpacking-arguments
        """

        print(self.vcf.header, end="")
        #print(*map(self.fix_record, self.vcf), sep="", end="")
        for line in self.vcf:
            print(self.fix_record(line), sep="", end="")

if __name__ == "__main__":
    app = Application()
    app.start()

Edited version of Namshin's code this can handle if REF or ALT is '-' ( these records in vcf file are observed when there is insertion/deletion).
Added documentation in the script to explain code block I have modified.
Hope this useful.
Tested example with GRCh37
python fix_vcf.py -f genome.fa -o out_fixed.vcf input.vcf

example input.vcf:
##fileformat=VCFv4.2
#CHROM  POS ID  REF ALT QUAL  FILTER  INFO
7 6027157 . - T 1 PASS  NS=1  GQ  10
9 98218616  . - C 1 PASS  NS=1  GQ  10
17  41245200  . AT  - 1 PASS  NS=1  GQ  10
13  32929213  . C - 1 PASS  NS=1  GQ  10

out_fixed.vcf:
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
#CHROM  POS ID  REF ALT QUAL  FILTER  INFO
7 6027157 . T TT  1 PASS  NS=1
9 98218616  . A AC  1 PASS  NS=1
17  41245199  . GAT G 1 PASS  NS=1
13  32929212  . CC  C 1 PASS  NS=1

# -*- coding: utf-8 -*-
import argparse
import sys

import pysam

class Application:
    def __init__(self):
        self.args = Application.get_args()
        self.vcf = pysam.VariantFile(sys.stdin)
        self.fasta = pysam.FastaFile(self.args.fasta)

    @staticmethod
    def get_args() -> argparse.Namespace:
        """
        Read command line arguments and set input and output sources.

        :return: parsed command line arguments
        """

        parser = argparse.ArgumentParser(prog="fixvcf.py", description="Fixes deletion represented by -")
        parser.add_argument('--version', action='version', version='%(prog)s 0.1')

        parser.add_argument("input", help="input file, use - to read from stdin")
        parser.add_argument("-f", "--fasta", help="fasta reference file", required=True)
        parser.add_argument("-o", "--output", help="output file")

        if len(sys.argv) == 1:
            parser.print_help(sys.stderr)
            sys.exit(1)

        args = parser.parse_args()

        if args.output:
            sys.stdout = open(args.output, "w")

        if args.input != "-":
            sys.stdin = open(args.input, "r")

        return args

    def fix_record(self, record: pysam.libcbcf.VariantRecord) ->     pysam.libcbcf.VariantRecord:
        """
        If the ALT column of the records contains a `-`, the variant position is shift one position to the left.
        For reaching this the reference base for the position is determined and prepended to the current REF bases.
        The ALT value `-` is replaced with the determined ref base on position before. In case there are multiple ALT
        values the determined ref base is prepended to them.

        If the REF column of the records contains a `-`, the variant position will change.
        For reaching this the reference base for the position is determined and prepended to the current ALT bases.
        The REF value `-` is replaced with the determined ref base on position. In case there are multiple ALT
        values the determined ref base is prepended to them.

        :param record: vcf file record
        :return: fixed vcf record
        """

        if "-" in record.alts:
            ref = self.fasta.fetch(reference=record.chrom, start=record.pos - 2, end=record.pos -1)

            record.pos = record.pos - 1
            record.ref = f"{ref}{record.ref}"
            record.alts = tuple(map(lambda alt: ref if alt == "-" else f"{ref}{alt}", record.alts))

        if "-" in record.ref:
            ref = self.fasta.fetch(reference=record.chrom, start=record.pos - 1, end=record.pos)

            record.pos = record.pos
            record.ref = f"{ref}"
            record.alts = tuple(map(lambda alt: f"{ref}{alt}", record.alts))

        return record

    def start(self):
        """
        Print the header of the vcf file. Afterwards the fix_record method is applied to each record in the vcf file.

        For more information about python's map() see:
https://docs.python.org/3.6/library/functions.html#map

        For more information about the unpacking mechanism using `*`see:
        https://docs.python.org/3/tutorial/controlflow.html#tut-unpacking-arguments
        """

        print(self.vcf.header, end="")
        #print(*map(self.fix_record, self.vcf), sep="", end="")
        for line in self.vcf:
            print(self.fix_record(line), sep="", end="")

if __name__ == "__main__":
    app = Application()
    app.start()


text:
After running a DRIMSeq pipeline, how do I know which genes are upregulated in the different conditions?
After running a DRIMSeq pipeline and obtaining the genes that are differently used between the null and full models, how do I select the genes that are differently used in the different conditions?
According to the vignette, I would have to use the plotProportions function and say in the group_variable the variable of my interest. What if I want a list of the genes that are differently used specifically in the variables of my interest? Checking every gene manually sounds exhausting. Help would be much appreciated.

After running the test for differential transcript usage you get a list with gene- and transcript-level p-values that you can use to filter which differences are significant. The gene-level p-value tells you if there is differential usage of transcripts within a gene and the transcript-level p-value tells you wich transcripts actually are changing.  Therefor you will want to select significant genes for further investigation.
The plots are mainly used to visualize top hits as detailed in the vignette.
res <- results(d)
# The data is ordered by p-value
res <- res[order(res$pvalue, decreasing = FALSE), ]
# The top hit is selected
top_gene_id <- res$gene_id[1]
plotProportions(d, gene_id = top_gene_id, group_variable = "group")

However you should consider running the two-stage test that is described on page 16 of the vignette you are using to reduce the number of false positive hits. From this you will get another list of adjusted p-values that can be used as above.

I'm very gratefull for your input. But what if I have more than one group_variable how do I know why p-values are related to that group_variable?

When you run the test (d <- dmTest(d, coef = "groupKD", verbose = 1)) you are specifying the contrast you are investigating (in this case "groupKD"). If you for example have two variables that you want to look at independently you will have to make two DRIMSeq data sets and run dMTest two times. It would be best if you can specify your model.matrix in the question so I can give you more specific answers.

Let's say I have a design matrix which is ~ treatment + celltype and my null model is ~1 . After running the pipeline with this design, how do I get a list of the genes whose transcripts are differentially used for the treatment or for the celltype ? I don't want to use two different models to get the DTU for each variable as the above additive model gave better gene precision estimates than models with individual variable.

If you have two columns in your design matrix you can use the coef approach. for example dmTest(d, coef = "treatmentB, verbose = 1) will give you the p values for comparing treatment B vs A (in this example A is the reference level, if you want to change that you can refactor the columns in the column data). But please note that in this case the null model becomes ~ cellltype which is what you need to isolate the treatment from the celltype effect.


text:
Argument stopped working
I am using ComplexHeatmap package for transcriptomic analyses and all of a sudden,  row_labels and column_labels as well as heatmap_width and heatmap_height stopped being recognized as arguments of the function Heatmap:
actual Heatmap
Heatmap(matrix [RowFunctionOrder,ColumnImposedOrder],
        col = colgrad,
        show_column_names = F,
        cluster_columns = FALSE,
        cluster_rows = FALSE,
        row_names_gp = gpar(fontsize=7),
        heatmap_width=1,
        heatmap_height = 1,
        row_labels =  rnaseqMatrix$Protein.ID,
        bottom_annotation = Column_ha
        )
#creating a column annotation

Column_ha = HeatmapAnnotation(condition = colnames(matrix),condition = runif(30), col = list(condition=Colcolour), show_legend = FALSE).

*Error in Heatmap(matrix[RowFunctionOrder, ColumnImposedOrder], col = colgrad,  :
  unused arguments (heatmap_width = 1, heatmap_height = 1, row_labels = rnaseqMatrix$Protein.ID)*.

this code was perfectly working this morning, I do not understand why suddenly the arguments were not valid anymore. they also seem to have disappeared from the help menu. this event followed a bug where I was receiving messaged saying the help(ComplexHeatmap) could't be update. this messages started appearing when I started typing arguments.

Hello Jordan, my best guess is that you have loaded a different package that has masked the original Heatmap function from ComplexHeatmap package. Does you have the same problem when you start a new R session and load only ComplexHeatmap package?

Please edit your question and add the output of sessionInfo(). Also, try using ComplexHeatmap::Heatmap() instead of just Heatmap() - that should tell you if you're running into what @KamilSJaron is talking about.