2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Michael Aug 13 '19 at 8:08
  • $\begingroup$ 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? $\endgroup$ – terdon Aug 13 '19 at 11:23
  • $\begingroup$ 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? $\endgroup$ – terdon Aug 14 '19 at 18:30
2
$\begingroup$

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – terdon Aug 13 '19 at 12:51
  • $\begingroup$ 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(). $\endgroup$ – haci Aug 13 '19 at 17:41
  • $\begingroup$ Thank you. That's very helpful, have an upvote! $\endgroup$ – terdon Aug 13 '19 at 18:00
1
$\begingroup$

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ I just tried this and (assuming I'm using the right executables, I found 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. $\endgroup$ – terdon Aug 13 '19 at 13:15
  • $\begingroup$ 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. $\endgroup$ – Hadi Aug 13 '19 at 17:18
  • $\begingroup$ 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. $\endgroup$ – Alex Reynolds Aug 13 '19 at 20:45
  • $\begingroup$ 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. $\endgroup$ – terdon Aug 14 '19 at 8:00
1
$\begingroup$

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – terdon Aug 14 '19 at 18:26
  • 1
    $\begingroup$ 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. $\endgroup$ – mRotten Aug 14 '19 at 18:38
  • $\begingroup$ 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 :) $\endgroup$ – terdon Aug 14 '19 at 19:41
  • $\begingroup$ @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. $\endgroup$ – mRotten Aug 15 '19 at 5:10
0
$\begingroup$

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Hadi Aug 13 '19 at 17:14
  • $\begingroup$ @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. $\endgroup$ – terdon Aug 13 '19 at 17:37
  • $\begingroup$ I have edited the original post for accurate examples and yes, my original files use tab as delimiter but output uses single space $\endgroup$ – Hadi Aug 14 '19 at 3:36
  • $\begingroup$ @Hadi see updated answer. It should work as you want now. $\endgroup$ – terdon Aug 14 '19 at 18:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.