Thursday, September 25, 2014

Autosomal Compare (*nix)

For comparing two autosomal files, there is a Windows tool called, Autosomal Segment Analyzer. However, in Unix-based systems, you can compare using the below command. The below command assumes that the autosomal files are in 23andMe format. If you have autosomal files in any other format, you should be able to convert using the commands provided in Autosomal Converter for *nix page

Prerequisites: Any Unix-based system

$ join --nocheck-order -e EMPTY --header file1.txt file2.txt |awk 'BEGIN { snp_threshold = 700; mb_threshold = 7; error_radius = 350; largest_mb = 0; total_mb = 0; count=0; seg_start=0;seg_end=0; chr=0; pchr=0; error_pos=0; print "\nChr\tStart Position\tEnd Position\tLen(Mb)\tSNPs";} { chr = $2; seg_len = (seg_end-seg_start)/1000000; if( !($4 == $7 || substr($4,1,1) == substr($7,1,1)|| substr($4,2,1) == substr($7,2,1) || substr($4,1,1) == substr($7,2,1)|| substr($4,2,1) == substr($7,1,1) ) || pchr!=chr) { if( seg_end - error_pos > error_radius ) { count++; seg_end = $3; } else { if( count > snp_threshold && seg_len > mb_threshold) { total_mb = total_mb + seg_len; if(largest_mb < seg_len) largest_mb = seg_len; print chr"\t"seg_start"\t"seg_end"\t"seg_len"\t"count; } count = 0; seg_start = $3; } error_pos=$3; } else { count++; seg_end = $3; } pchr = chr;}END {print "\nLargest Segment: "largest_mb" Mb";print "Total Shared: "total_mb" Mb\n";}'

Note: file1.txt and file2.txt are the two files being compared. The SNP threshold of 700, Mb Threshold of 7 and error radius of 350 SNPs can be modified which are bolded for convenience.

Screenshot: