projects/radseq/cmd.lst

   1 rm -fr work
   2 mkdir lst work
   3 #find ./fq.p -name '*.gz'|sort > ./lst/rawfqp.lst
   4 #find ./fq.f -name '*.gz'|sort > ./lst/rawfqf.lst
   5 find /share/users/huxs/work/tiger/fq.p /share/users/huxs/work/tiger/fq.f -name '*.gz'|sort > ./lst/rawfq0.lst
   6 mkdir work/radseq work/parents
   7 ./src/rawfq2list.pl ./lst/rawfq0.lst ./lst/rawfq.lst
   8
   9 make -f make/makefile
  10 #ps -ef|grep cutadapt|sort -rnk4|les -N
  11
  12 ====================================
  13
  14 ./src/separatebytag.pl tigrad.lst 0,0,0,0 work/radseq/A1_1.1.cut.gz work/radseq/A1_1.2.cut.gz work/f1/A1_1_0000
  15 ./src/separatebytag.pl tigrad.lst 0,0,0,0 work/radseq/B1_2.1.cut.gz work/radseq/B1_2.2.cut.gz work/f1/B1_2_0000 &
  16 ./src/separatebytag.pl tigrad.lst 0,0,0,0 work/radseq/C1_3.1.cut.gz work/radseq/C1_3.2.cut.gz work/f1/C1_3_0000 &
  17 ./src/separatebytag.pl tigrad.lst 0,0,0,0 work/radseq/d1_4.1.cut.gz work/radseq/d1_4.2.cut.gz work/f1/d1_4_0000 &
  18
  19 ./src/separatebytag.pl tigrad.lst 210210 work/f1/A1_1_0000.NA.Unknown.1.fq.gz work/f1/A1_1_0000.NA.Unknown.2.fq.gz work/f1/A1_1_210210 &
  20 ./src/separatebytag.pl tigrad.lst 210210 work/f1/B1_2_0000.NA.Unknown.1.fq.gz work/f1/B1_2_0000.NA.Unknown.2.fq.gz work/f1/B1_2_210210 &
  21 ./src/separatebytag.pl tigrad.lst 210210 work/f1/C1_3_0000.NA.Unknown.1.fq.gz work/f1/C1_3_0000.NA.Unknown.2.fq.gz work/f1/C1_3_210210 &
  22 ./src/separatebytag.pl tigrad.lst 210210 work/f1/d1_4_0000.NA.Unknown.1.fq.gz work/f1/d1_4_0000.NA.Unknown.2.fq.gz work/f1/d1_4_210210 &
  23
  24 ./src/separatebytag.pl tigrad.lst 221210 work/f1/A1_1_210210.NA.Unknown.1.fq.gz work/f1/A1_1_210210.NA.Unknown.2.fq.gz work/f1/A1_1_221210 &
  25 ./src/separatebytag.pl tigrad.lst 221210 work/f1/B1_2_210210.NA.Unknown.1.fq.gz work/f1/B1_2_210210.NA.Unknown.2.fq.gz work/f1/B1_2_221210 &
  26 ./src/separatebytag.pl tigrad.lst 221210 work/f1/C1_3_210210.NA.Unknown.1.fq.gz work/f1/C1_3_210210.NA.Unknown.2.fq.gz work/f1/C1_3_221210 &
  27 ./src/separatebytag.pl tigrad.lst 221210 work/f1/d1_4_210210.NA.Unknown.1.fq.gz work/f1/d1_4_210210.NA.Unknown.2.fq.gz work/f1/d1_4_221210 &
  28
  29
  30 Parents WGS:
  31 bwa aln -l 17 -q 10
  32
  33 RADseq:
  34 bwa aln -l 17 -B 5
  35 【only read1. read2 still with -q 10】
  36
  37 __BIOPIC__
  38 cd /home/huxuesong/analysis/tiger/20120712
  39 find ./fq/parents/ -name '*.1.cut.gz'|perl -lane 's/\.1\.cut\.gz$//;print' > parents.lst
  40 ls -1 ./fq/f1/*.1.fq.gz|perl -lane 's/\.1\.fq\.gz$//;print'|sort > f1.lst
  41
  42 soapsnp -u -t -L 101 -z '!' -Q 41 -d ref/tigris.fa
  43
  44 ./fq/parents/BHX011_11
  45 ./fq/parents/BHX011_11.1.cut.gz
  46 ./fq/parents/BHX011_11.2.cut.gz
  47
  48
  49 __SH__
  50 #!/bin/bash
  51 #$ -v PERL5LIB,PATH,PYTHONPATH,LD_LIBRARY_PATH
  52 #$ -cwd -r y -l vf=7.5g,p=2
  53 #$ -o /dev/null -e /dev/null
  54 #$ -S /bin/bash -t 1-2
  55
  56 SEEDFILE=parents.lst
  57 SEED=$(sed -n -e "${SGE_TASK_ID} p" $SEEDFILE)
  58 MAIN=`basename "$SEED"`
  59
  60 bwa aln -Y -l 17 -q 10 -t 6 tigris ${SEED}.1.cut.gz >./sai/${PPID}_${JOB_ID}_${MAIN}_1.sai 2>./sai/${PPID}_${JOB_ID}_${MAIN}_1.log
  61 bwa aln -Y -l 17 -q 10 -t 6 tigris ${SEED}.2.cut.gz >./sai/${PPID}_${JOB_ID}_${MAIN}_2.sai 2>./sai/${PPID}_${JOB_ID}_${MAIN}_2.log
  62
  63 bwa sampe -a 800 tigris ./sai/${PPID}_${JOB_ID}_${MAIN}_1.sai ./sai/${PPID}_${JOB_ID}_${MAIN}_2.sai ${SEED}.1.cut.gz ${SEED}.2.cut.gz 2>./bwatg/${MAIN}.samlog | gzip -9c >./bwatg/${MAIN}.sam.gz
  64 date >> ./bwatg/${MAIN}.samlog
  65 echo "sai: ${PPID}_${JOB_ID}_${MAIN}" >> ./bwatg/${MAIN}.samlog
  66
  67 samtools view -bS ./bwatg/${MAIN}.sam.gz > ./bwatg/${MAIN}.bam 2>./bwatg/${MAIN}.bamlog
  68
  69 samtools sort ./bwatg/${MAIN}.bam ./bwatg/${MAIN}.sort
  70
  71 samtools index ./bwatg/${MAIN}.sort.bam &
  72
  73
  74 [after bcf]
  75 ./bcf2ped.pl radseq.tfam radseq.bcgv.bcf radseq17
  76
  77 p-link --tfile radseq17 --reference-allele radseq17.MinorAllele --fisher --out Pradseq17 --model --cell 0
  78
  79 grep REC Pradseq17.model > Pradseq17.model.REC.0
  80
  81 # p-link output may lack A2
  82 perl -lane 'if (@F<8) {splice @F,3,0,".";} die if @F<8; print join("\t",@F)' Pradseq17.model.REC.0 > Pradseq17.model.REC
  83
  84 sort -k 8,8n -k 2,2 Pradseq17.model.REC > Pradseq17.model.REC.nk8
  85 ./dojoin.pl radseq17.dict 3 Pradseq17.model.REC.nk8 2 tmp.REC
  86
  87 #awk '{print $3" "$1" "$2" "$4" "$6" "$7" "$8" "$9" "$10" "$11}' tmp.REC.out > tmp.REC.out.j
  88 awk 'BEGIN {OFS="\t"} {print $3,$1,$2,$4,$6,$7,$8,$9,$10,$11}' tmp.REC.out > tmp.REC.out.j
  89
  90 sort -k 10,10n -k 1,1 tmp.REC.out.j > tmp.REC.out.s
  91
  92 cat tmp.REC.out.j|perl -lane '@a=split /\//,"$F[7]/$F[8]";$sum=$a[0]+$a[1]+$a[2]+$a[3];
  93 $theta=($a[1]+$a[2])/($sum*2);
  94 $p=((1-$theta)**(2*$sum-($a[1]+$a[2])))*($theta**($a[1]+$a[2]))/(0.5**(2*$sum));
  95 $LOD=int(0.5+log($p)*1000/log(10))/1000;
  96 print join("\t",@F,$sum,int(0.5+$theta*100000)/100000,int($LOD),$LOD)' |sort -k2,2 -k3,3n > rec.npa
  97
  98 sort -k13,13n -k2,2 -k3,3n  rec.npa > rec.npas
  99
 100 ./plotscaf.pl rec.npa tig2cat
 101 ./plotscaf.pl rec.npa tig2cats
 102
 103 ./ped2hap.pl tighap.lst radseq17 dohap17
 104
 105 grep -P 'scaffold(75|1458|188)\t' rec.npa > rec.npa.sss
 106
 107 perl -lane 'print $_ if $F[-4]>15' rec16q20.npa.sss > rec16q20.npa.sss16
 108
 109 __END__
 110
 111 $ git push
 112 Counting objects: 11, done.
 113 Delta compression using up to 24 threads.
 114 Compressing objects: 100% (6/6), done.
 115 Writing objects: 100% (6/6), 938 bytes, done.
 116 Total 6 (delta 4), reused 2 (delta 0)
 117 error: file write error (No space left on device)
 118 fatal: unable to write sha1 file
 119 error: unpack failed: unpacker exited with error code
 120 To ssh://galaxy@repo.or.cz/srv/git/GalaxyCodeBases.git
 121  ! [remote rejected] master -> master (n/a (unpacker error))
 122 error: failed to push some refs to 'ssh://galaxy@repo.or.cz/srv/git/GalaxyCodeBases.git'
 123
 124 这是远端故障……
 125 @ 10:23:26 Tue Sep 11 2012