ASE分析
1、Prepare necessary input files（可以参考上次的博客http://www.cnblogs.com/renping/p/7391028.html）
1）对fq1和fq2合并
cat fq1 fq2
2）对bam 文件转换成psl格式
/share/nas2/genome/biosoft/Python/2.7.8/bin/python /share/nas1/wenyh/develop/tools/Au-public-master/iron/utilities/sam_to_psl.py -r transcript.fa T16.bam >T16.psl
3）gtf format convert to gpd format
/share/nas1/wenyh/develop/tools/gtfToGenePred transcript.gtf -genePredExt transcript.gpd.tmp
awk '{print 0"\t"$0}'transcript.gpd.tmp >transcript.gpd.tmp2
/share/nas1/wenyh/develop/pacbio/IDP-ASE/julia/bin/julia /home/wenyh/.julia/v0.4/IDPASE/scripts/convert_gpd.jl transcript.gpd.tmp2 >transcript.gpd.tmp3
4）vcf注释和选杂合的vcf文件
注释vcf文件。（参考博客：http://www.cnblogs.com/renping/p/7467348.html）
awk '$10!~/1\/1/;$10!~/\.\/\./{print}'|le >final.snp.anno.vcf1 ##筛选杂合
le final.snp.anno.vcf1|grep -v '#'|cut -f 1 |sort |uniq -c | awk '{print $2,$1}'|less -S|sort -k 2nr|le >Snp.distribution
2、Prepare Gene level data
1) mkdir temp/; mkdir gene_files; mkdir isoform_files; mkdir gene_out; mkdir isoform_out;
2) for i in `le snp.distribution |awk '$1<10 {print $2}'|le`; do echo "/share/nas1/yangch/tools/julia/bin/julia -p 4 /home/yangch/.julia/v0.4/IDPASE/src/prep_runs.jl \
-a /share/nas1/yangch/RENPP/out/T19.psl \
-g /share/nas1/yangch/RENPP/out/transcript.gpd.tmp3 \
-v /share/nas1/yangch/RENPP/out/final.snp.anno.vcf1 \
-q /share/nas1/yangch/RENPP/out/T19.fq \
-d /share/nas1/yangch/RENPP/out/temp \
-c ${i} \
-f 1 \
-o /share/nas1/yangch/RENPP/out/gene_files/ \
-p T19 "; done >A1.sh #####Prepare Gene level data
3) for i in `ls /share/nas1/yangch/RENPP/out/gene_files/|perl -lne '{next if /^\s+$/;/T19_(reads|true)_(.*)\.txt/;print $2}'|sort|uniq|less`;\
do echo "/share/nas1/yangch/tools/julia/bin/julia -p 4 /home/yangch/.julia/v0.4/IDPASE/scripts//phase_by_loci_sub.jl \
-t /share/nas1/yangch/RENPP/out/gene_files/T19_true_${i}.txt \
-a /share/nas1/yangch/RENPP/out/gene_files/T19_reads_${i}.txt \
-o /share/nas1/yangch/RENPP/out/gene_out/ \
-l 1 \
-r ${i} \
-i 10000 \
-b 1000 \
-c 4 \
-d /home/yangch/.julia/v0.4/IDPASE/scripts/ \
-n SGS \
-m 1 0 \
-s 1.0"; done >to_run_curr.sh #### Get commands to run each gene individually
4) Concatenate all gene level results
find gene_out/ -name "REAL*" | xargs cat > gene_out/gene_results.txt
