Merging two .txt files

Question

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.

Answer 1

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.

Answer 2

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.

Answer 3

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.

Answer 4

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.