reformat_cnv.sh

for chr in 10 # X # 2 6 11 # {7..20} # {2..20} # {1..20} #22 # {1..8} # {1..22} X MT XY Y
do
echo $chr
date
 mkdir temporary_files/chr"$chr"_baf1
 mkdir temporary_files/chr"$chr"_baf2
 mkdir temporary_files/chr"$chr"_baf3
 mkdir temporary_files/chr"$chr"_baf4
 mkdir temporary_files/chr"$chr"_l2r1
 mkdir temporary_files/chr"$chr"_l2r2
 mkdir temporary_files/chr"$chr"_l2r3
 mkdir temporary_files/chr"$chr"_l2r4
 mkdir chr"$chr"_complete
  
##Split the file by 100 lines at a time
#l2r
 
 echo "splitting l2r"
 date
 zcat ukb_l2r_chr"$chr"_v2.txt.gz | split -l 100 --suffix-length=6  -d - temporary_files/chr"$chr"_l2r1/chr"$chr"_l2r_
#baf
 echo "splitting baf"
 date
 zcat ukb_baf_chr"$chr"_v2.txt.gz | split -l 100 --suffix-length=6  -d - temporary_files/chr"$chr"_baf1/chr"$chr"_baf_

 #Benchmark: with 100 per file, can do about 400 SNPs/minute. For a whole genome thats ~1500 minutes (roughly 24 hours). So just over 2 days to do this for all sets
 #With 50 per file, same speed. No speed up here... 
 #With 25 per file, same speed
 #with 5 per file, 73 files per minute = 365 snps/minute. i dont think this isa  slow down, just timing is more accurate given file inputes
 #with 1 per file, 359 per minute.. so no real timing change here (6 SNPs per second).
 
##Transpose (I think this is the long step)

#l2r 
 filecount=1
 for files in $(ls temporary_files/chr"$chr"_l2r1)
  do
  perl transpose_space.pl  temporary_files/chr"$chr"_l2r1/$files > temporary_files/chr"$chr"_l2r2/$files 
  if [ $(($filecount % 100))  -eq 0 ]
  then 
   echo "l2r is at $filecount"
   date
  fi
   filecount=$(($filecount+1))
 done
 
#baf

filecount=1
 for files in $(ls temporary_files/chr"$chr"_baf1)
 do
  perl transpose_space.pl  temporary_files/chr"$chr"_baf1/$files > temporary_files/chr"$chr"_baf2/$files 

  if [ $(($filecount % 100))  -eq 0 ]
  then 
   echo "baf is at $filecount"
   date
  fi
   filecount=$(($filecount+1))
 done

 #with 100 lines, 1 transpose took 50 seconds. So again, 2 days to do all of this
 #with 50 lines, 1 transpose took 20 seconds. 2 took 40 seconds.. so potential speed gain here. (150/min)
 #With 25 lines, 1 transpose took..9 seconds..so the speedup isnt great here.. (150/min)
 #with 5 lines, 2.2 seconds... 25 lines might be best. 
 #With a sed based transpose, 230 files in 1 minute...so this is much faster than the other approaches, provided I can paste themm....another says 85 in 30 seconds. so thats 170/min.. really not much faster?? (170/min)

##Rejoin file 

#L2r
echo "Rejoining l2r"
date
time paste temporary_files/chr"$chr"_l2r2/* > temporary_files/chr"$chr"_l2r3/ukb_l2r_chr"$chr"_v2_t.txt
#baf
echo "Rejoining baf"
date
paste temporary_files/chr"$chr"_baf2/* > temporary_files/chr"$chr"_baf3/ukb_baf_chr"$chr"_v2_t.txt

#WAs 62 minutes with 50 lines per file.

#Split file by line (I could alternatively paste after this but its too many i/o
#l2r
echo "spliting l2r"
date
split -l 1  --suffix-length=7  -d temporary_files/chr"$chr"_l2r3/ukb_l2r_chr"$chr"_v2_t.txt temporary_files/chr"$chr"_l2r4/ukb_l2r_chr"$chr"_
#baf
echo "spliting baf"
date
split -l 1  --suffix-length=7  -d temporary_files/chr"$chr"_baf3/ukb_baf_chr"$chr"_v2_t.txt temporary_files/chr"$chr"_baf4/ukb_baf_chr"$chr"_

Based on previous stuff, this will take 2 days
rm -rf temporary_files/chr"$chr"_l2r1
rm -rf temporary_files/chr"$chr"_baf1
rm -rf temporary_files/chr"$chr"_l2r2
rm -rf temporary_files/chr"$chr"_baf2
rm -rf temporary_files/chr"$chr"_l2r3
rm -rf temporary_files/chr"$chr"_baf3

#For each subject:
 #perl transpose_space.pl <( 
 #Concatenate together genotype Names  and transpose back..
 echo "adjusting bim file"
 date
 perl transpose.pl <(cat ukb_snp_chr"$chr"_v2.bim | awk 'BEGIN{OFS="\t"}{print $2}')  > chr"$chr"_pcv_snps.txt
 cat ukb_snp_chr"$chr"_v2.bim | awk 'BEGIN{OFS="\t"}{print $2}'  > chr"$chr"_pcv_snps_ALT.txt
 

done


#need this for chr 2, 6?

for chr in {1..22} #X x
do
 awk 'BEGIN{OFS="\t"}{print $2,$1,$4}'  ukb_snp_chr"$chr"_v2.bim  > ukb_snp_chr"$chr"_v2.pcv_positions
done
echo -e "SNP Name\tChr\tPosition" > posheader.txt
 cat posheader.txt  *.pcv_positions > allchr.pcvpos

#Get priority list of IDs
seq -w 0 488377   
LC_ALL=C join <(LC_ALL=C sort -k1b,1 UKB_ptsd_eur_related_m3_may3_2019.pheno | awk '{print $1}' ) <( paste <(awk '{print $1}' /mnt/ukbb/adam/ptsd/preimputation_genotypes/ukb41209_cal_chr1_v2_s488292.fam) <(seq -w 0 488377 ) | LC_ALL=C sort -k1b,1) | sort -k2 > priority_subjects.txt
LC_ALL=C join <(LC_ALL=C sort -k1b,1 UKB_ptsd_eur_related_m3_may3_2019.pheno | awk '{print $1}' ) <( paste <(awk '{print $1}' /mnt/ukbb/adam/ptsd/preimputation_genotypes/ukb41209_cal_chr1_v2_s488292.fam) <(seq -w 0 488377 ) | LC_ALL=C sort -k1b,1) | sort -k2 > priority_subjects.txt




split --numeric-suffixes=1 -l 5000 priority_subjects.txt priority
 
 mkdir subjects_temporary
 mkdir subjects_complete
 
    # if [ ! -f "chr"$chr"_complete/"$chr"_"$subject".pcv" ]
   # then
   # fi
     # mv temporary_files/chr"$chr"_l2r4/ukb_l2r_chr"$chr"_0$subject temporary_files/chr"$chr"_l2r4_"$n"/.
   #mv temporary_files/chr"$chr"_baf4/ukb_baf_chr"$chr"_0$subject temporary_files/chr"$chr"_baf4_"$n"/.
   
   #check all chromosomes now...
   for chr in {1..22}
   do
    
    ls --sort=non temporary_files/chr"$chr"_baf4 | wc -l 
     ls --sort=non temporary_files/chr"$chr"_l2r4 | wc -l 
    
    done
    
    bafs fromchromsome 10 are missing 123 subjects... fuck them...
    
for n in   26 27 28 29 30 # 01 02 03 04 05 06 07 08 09 10 11 12  13 14 15 16 # 17 18 19 20 21 22 23 24 25
do
 now=$(date)
 echo "It is $now. Processing chunk $n"
 
 for subject in 080000 #$(awk '{print $2}' priority$n )  # go to number of subjects 488377
 do
  if [ ! -f "subjects_complete/"$subject".pcv" ]
  then
  
   for chr in  {1..22} #X #
   do
     paste chr"$chr"_pcv_snps_ALT.txt <(sed 's/\t/\n/g' temporary_files/chr"$chr"_baf4/ukb_baf_chr"$chr"_0$subject)  <(sed 's/\t/\n/g' temporary_files/chr"$chr"_l2r4/ukb_l2r_chr"$chr"_0$subject)  >  subjects_temporary/"$chr"_"$subject".pcv 
   done
   
   echo -e "Name\t$subject.B Allele Freq\t$subject.Log R Ratio" > subjects_temporary/00_"$subject".pcv
   cat subjects_temporary/*_"$subject".pcv  | grep -v NA > subjects_complete/"$subject".pcv
   rm -f  subjects_temporary/*_"$subject".pcv
  fi
 done
 echo "It is $now. Finished chunk $n"
done

#figureo ut what happened tp /oasis/scratch/comet/amaihofer/temp_project/plate29/453044.pcv: line 243389  
subject=453044
wc -l subjects_temporary/*_"$subject".pcv #problem happens at line 243.. so just fine out where cat puts that...
#I think chromosmoe 1... ok... fucker... asshole
... yep its fucked somehow... too many lines compared to rsids...
#look at the tail of each file...
tail -n1 subjects_temporary/*_453044.pcv
#Ok the problem lies with just chromosome 1... depressing!
#this has 95224 in the file, in the bim there are 63K. what the fuckity fuck. Check the unmodified original file.


#note some parts of chr 10 and 11 have been moved into folders temporary_files/chr"$chr"_l2r4_"$n"/. temporary_files/chr"$chr"_baf4_"$n"/. (possibly also 9) do not retain this systems
#chr 13 and 14 have done 1:10 so far
#5 chr15 has done5-10

##Fingerprinting 
awk 'NR>=450000&&NR<=450200{print $1}' /mnt/ukbb_cnv/CNV_example/subjects_complete/000000.pcv  > fingerprint.snplist
awk 'NR>=450000&&NR<=450200' /mnt/ukbb_cnv/CNV_example/subjects_complete/269225.pcv > fingerprint.penn

awk 'NR==269226{print $1,$2}' /mnt/ukbb/ancestry_unzip/ukb_chr2_v2.fam > fingerprint_test.subject #subject counting starts at 0 for fingerprint files, need subject minus 1
plink --bfile /mnt/ukbb/ancestry_unzip/ukb_chr10_v2 --keep fingerprint_test.subject --extract fingerprint.snplist --recodeA --out fingerprint.alleles


/mnt/ukbb_cnv/CNV_example/temporary_files/chr10_baf1/

#nEXT: Then concatenate all files.. and add a header

#Get first 100 SNps for first subject.. ok this i know fingerprints correctly, which is good..tail 
 awk '{print $1}' /mnt/ukbb_cnv/CNV_example/temporary_files/chr10_baf1/chr10_baf_000000  > baf.snp
 head -n100 /mnt/ukbb/ancestry_unzip/ukb_chr10_v2.bim > first100snps.snplist
 awk 'NR==1 {print $1,$2}' /mnt/ukbb/ancestry_unzip/ukb_chr1_v2.fam > fingerprint_test.subject 
 plink --bfile /mnt/ukbb/ancestry_unzip/ukb_chr10_v2 --keep fingerprint_test.subject --extract first100snps.snplist --recodeA --out fingerprint.alleles

 #compare also to the finished file
 head /mnt/ukbb_cnv/CNV_example/subjects_complete/000000.pcv 
awk 'NR>=450000&&NR<=450200{print $1}' /mnt/ukbb_cnv/CNV_example/subjects_complete/000000.pcv  > fingerprint.snplist

#works on random people... check some established files to be sure, but i think this is it
tail -n500 /mnt/ukbb_cnv/CNV_example/subjects_complete/419835.pcv > baf.comp
tail -n500  /mnt/ukbb_cnv/CNV_example/subjects_complete/419835.pcv | awk '{print $1}'  > fingerprint2.snplist
awk 'NR==419836 {print $1,$2}' /mnt/ukbb/ancestry_unzip/ukb_chr1_v2.fam > fingerprint_test.subject  
 plink --bfile /mnt/ukbb/ancestry_unzip/ukb_chr9_v2 --keep fingerprint_test.subject --extract  fingerprint2.snplist --recodeA --out fingerprint.alleles2



awk 'NR>=350000&&NR<=350200{print }'  /mnt/ukbb_cnv/CNV_example/subjects_complete/419835.pcv > baf.comp
awk 'NR>=350000&&NR<=350200{print $1}' /mnt/ukbb_cnv/CNV_example/subjects_complete/419835.pcv  > fingerprint.snplist
awk 'NR==419836 {print $1,$2}' /mnt/ukbb/ancestry_unzip/ukb_chr1_v2.fam > fingerprint_test.subject  
 plink --bfile /mnt/ukbb/ancestry_unzip/ukb_chr19_v2 --keep fingerprint_test.subject --extract  fingerprint.snplist --recodeA --out fingerprint.alleles2

#damn im good..




1000 subjects in 1 minute for 12000 SNPs... 2 minutes for 100000 SNPs...so scaling is not linear...1 day to finish..


so 400 minutes per thing... thats 7 hours for just one chromosome... 20% faster (600 in 30 seconds) if i didnt have to fuck with the headers, which shouldnt be in at this point anyway because im going to
can do 2000 in 1 minute using the non tranpose way...




#join with SNP names as well as BAFs processed in the the same way

#file checking info

for chr in  {1..22}  # 21 22
do
 
# echo chr$chr baf1; 
#ls --sort=non temporary_files/chr"$chr"_baf1  | wc -l
# echo chr$chr baf2; 
#ls --sort=non temporary_files/chr"$chr"_baf2 | wc -l
# echo chr$chr  baf3; 
#ls --sort=non temporary_files/chr"$chr"_baf3 | wc -l
 echo chr$chr  baf4; 
ls --sort=non temporary_files/chr"$chr"_baf4 | wc -l


# echo chr$chr l2r1; 
#ls --sort=non temporary_files/chr"$chr"_l2r1 | wc -l
# echo chr$chr l2r1; 
#ls --sort=non temporary_files/chr"$chr"_l2r2 | wc -l
# echo chr$chr  l2r1; 
#ls --sort=non temporary_files/chr"$chr"_l2r3 | wc -l
echo chr$chr  l2r4; 
ls --sort=non temporary_files/chr"$chr"_l2r4 | wc -l

done


chr10 missing for the following subjkects
sed: can't read temporary_files/chr10_baf4/ukb_baf_chr10_0488362: No such file or directory
sed: can't read temporary_files/chr10_baf4/ukb_baf_chr10_0488363: No such file or directory
sed: can't read temporary_files/chr10_baf4/ukb_baf_chr10_0488369: No such file or directory
sed: can't read temporary_files/chr10_baf4/ukb_baf_chr10_0488370: No such file or directory
sed: can't read temporary_files/chr10_baf4/ukb_baf_chr10_0488372: No such file or directory
sed: can't read temporary_files/chr10_baf4/ukb_baf_chr10_0488373: No such file or directory
sed: can't read temporary_files/chr10_baf4/ukb_baf_chr10_0488374: No such file or directory
sed: can't read temporary_files/chr10_baf4/ukb_baf_chr10_0488375: No such file or directory