API tools faq
paste
Guest User

Untitled

a guest
Jul 16th, 2019
90
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 7.59 KB | None | 0 0
raw download clone embed print report
  1. #!/bin/sh
  2. prefix=$(echo "$1" | cut -f 1 -d '-')
  3. cores=$(nproc --all)
  4. if [ "$1" = "" ]
  5. then
  6. echo "Llamada incorrecta, debe indicar el nombre del archivo Input en la llamada al script."
  7. echo "Por ejemplo: ./script.sh 20190618-0-Input_R1_.fastq.gz"
  8. echo "Se debe poner los archivos de entrada en la ruta del usuario actual: ${HOME}"
  9. else
  10. if [ "$2" = "" ]
  11. then
  12. echo ""
  13. else
  14. cores=$2
  15. fi
  16. cd ${HOME}
  17. if [ $(dpkg-query -W -f='${Status}' samtools 2>/dev/null | grep -c "ok installed") -eq 0 ];
  18. then
  19. sudo apt install samtools;
  20. fi
  21. if [ $(dpkg-query -W -f='${Status}' tabix 2>/dev/null | grep -c "ok installed") -eq 0 ];
  22. then
  23. sudo apt install tabix;
  24. fi
  25. if [ $(dpkg-query -W -f='${Status}' bowtie2 2>/dev/null | grep -c "ok installed") -eq 0 ];
  26. then
  27. sudo apt install bowtie2;
  28. fi
  29. if [ $(dpkg-query -W -f='${Status}' curl 2>/dev/null | grep -c "ok installed") -eq 0 ];
  30. then
  31. sudo apt install curl;
  32. fi
  33. if [ $(dpkg-query -W -f='${Status}' bcftools 2>/dev/null | grep -c "ok installed") -eq 0 ];
  34. then
  35. sudo apt install bcftools;
  36. fi
  37.  
  38. if test -f "$prefix-0-Input_R1_.fastq.gz"; then
  39. if test -f "bin/trimmomatic.jar"; then
  40. echo ""
  41. else
  42. echo " === Error en el preprocesado. El archivo bin/trimmomatic.jar no existe. Descargando... ==="
  43. curl -L https://github.com/timflutre/trimmomatic/archive/master.zip -o master.zip
  44. unzip master.zip
  45. cd trimmomatic-master
  46. make
  47. make install
  48. cp adapters/TruSeq3-SE.fa ${HOME}/
  49. cp adapters/TruSeq3-PE.fa ${HOME}/
  50. echo "=== Instalacion de Trimmomatic finalizada. ==="
  51. cd ${HOME}
  52. fi
  53. echo "=== Realizando preprocesado ==="
  54. echo ""
  55. echo "Se usaran $cores threads. Para cambiar este numero introduzca el numero de threads en la ejecucion del script, por ejemplo ./script.sh 20190618-0-Input_R1_.fastq.gz 4"
  56. echo ""
  57. echo ""
  58.  
  59. java -jar bin/trimmomatic.jar PE -basein $prefix-0-Input_R1_.fastq.gz -baseout $prefix-1-TRIM.fastq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3
  60.  
  61.  
  62.  
  63. if test -f "hg19.fa"; then
  64. echo ""
  65. else
  66.  
  67. if test -f "hg19.fa.gz"; then
  68. echo "=== Se ha encontrado hg19.fa.gz, descomprimiendo... ==="
  69. gunzip hg19.fa.gz
  70. else
  71. echo ""
  72. echo "No se ha encontrado el archivo hg19.fa, descargando..."
  73. echo ""
  74. wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/hg19.fa.gz
  75. gunzip hg19.fa.gz
  76. fi
  77. fi
  78. echo ""
  79. echo "=== Realizando alineamiento ==="
  80. echo ""
  81. if test -f "index/hg19_index.1.bt2"; then
  82. echo ""
  83. else
  84. echo ""
  85. echo "=== Archivos de referencia no encontrados. Este paso puede tardar unas horas ==="
  86. echo ""
  87.  
  88. mkdir index
  89. bowtie2-build hg19.fa.gz index/hg19_index
  90. fi
  91.  
  92. bowtie2 -x index/hg19_index -p $cores -1 $prefix-1-TRIM_1P.fastq.gz -2 $prefix-1-TRIM_2P.fastq.gz -S $prefix-2-BOWTIE.sam
  93.  
  94. echo ""
  95. echo "=== Ordenando ==="
  96. echo ""
  97. samtools sort -@ $cores $prefix-2-BOWTIE.sam > $prefix-3-SAMTOOLSsort.bam
  98. echo ""
  99. echo "=== Indexando ==="
  100. echo ""
  101. samtools index -@ $cores $prefix-3-SAMTOOLSsort.bam
  102. echo "Indexado terminado."
  103. echo ""
  104. echo "=== Opcion BCFTOOLS: ==="
  105. echo ""
  106.  
  107. bcftools mpileup -f hg19.fa $prefix-3-SAMTOOLSsort.bam | bcftools call -mv -Oz -o $prefix-4a-BCF.vcf.gz
  108.  
  109. bcftools filter -threads $cores -i 'QUAL>=20 & DP>=10' $prefix-4a-BCF.vcf.gz > $prefix-4a-BCF_filt.vcf.gz
  110.  
  111. echo ""
  112. echo "=== Opcion VARSCAN ==="
  113. echo ""
  114.  
  115. if test -f "VarScan.v2.4.3.jar"; then
  116. echo ""
  117. else
  118. echo ""
  119. echo "=== No se ha encontrado el archivo VarScan.v2.4.3.jar, descargando... ==="
  120. echo ""
  121. wget https://github.com/dkoboldt/varscan/raw/master/VarScan.v2.4.3.jar
  122. fi
  123.  
  124. echo "SNP + INDEL Variants"
  125.  
  126. samtools mpileup -f hg19.fa $prefix-3-SAMTOOLSsort.bam | java -jar VarScan.v2.4.3.jar mpileup2cns -v --min-coverage 100 --min-reads2 4 --min-avg-qual 20 --min-var-freq 0.30 --output-vcf --variants > $prefix-4b-VSCANcns_filt.vcf
  127.  
  128. echo "SNP Variants"
  129. samtools mpileup -f hg19.fa $prefix-3-SAMTOOLSsort.bam | java -jar VarScan.v2.4.3.jar mpileup2snp -v --min-coverage 100 --min-reads2 4 --min-avg-qual 20 --min-var-freq 0.30 --output-vcf --variants > $prefix-4b-VSCANsnp_filt.vcf
  130. echo "INDEL Variants"
  131. samtools mpileup -f hg19.fa $prefix-3-SAMTOOLSsort.bam | java -jar VarScan.v2.4.3.jar mpileup2indel -v --min-coverage 100 --min-reads2 4 --min-avg-qual 20 --min-var-freq 0.30 --output-vcf --variants > $prefix-4b-VSCANindel_filt.vcf
  132.  
  133. echo ""
  134. echo "=== Replacement of read groups in BAM file ==="
  135. echo ""
  136. if test -f "gatk-4.1.2.0/gatk-package-4.1.2.0-local.jar"; then
  137. echo ""
  138. else
  139. echo ""
  140. echo "=== No se ha encontrado el archivo gatk-package-4.1.2.0-local.jar, descargando... ==="
  141. echo ""
  142. wget https://github.com/broadinstitute/gatk/releases/download/4.1.2.0/gatk-4.1.2.0.zip
  143. unzip gatk-4.1.2.0.zip
  144. export PATH=$PATH:${HOME}/gatk-4.1.2.0/
  145. echo ""
  146. echo "Instalacion de GATK finalizada..."
  147. echo ""
  148. fi
  149.  
  150.  
  151. if test -f "All_20180423.vcf.gz"; then
  152. echo ""
  153. else
  154. echo ""
  155. echo "=== Archivo All_20180423.vcf.gz no encontrado, descargando... ==="
  156. echo ""
  157.  
  158. wget ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh37p13/VCF/All_20180423.vcf.gz
  159.  
  160. fi
  161.  
  162.  
  163. if test -f "All_20180423.vcf.gz.tbi"; then
  164. echo ""
  165. else
  166. echo ""
  167. echo "=== Archivo All_20180423.vcf.gz.tbi no encontrado, descargando... ==="
  168. echo ""
  169.  
  170. tabix -p vcf All_20180423.vcf.gz
  171.  
  172. fi
  173.  
  174. java -jar gatk-4.1.2.0/gatk-package-4.1.2.0-local.jar AddOrReplaceReadGroups -I $prefix-3-SAMTOOLSsort.bam -O $prefix-4c-GATKgroup.bam -ID 1 -LB dna -PL illumina -PU UNKNOWN -SM $prefix-3-SAMTOOLSsort
  175.  
  176. java -jar gatk-4.1.2.0/gatk-package-4.1.2.0-local.jar CreateSequenceDictionary -R hg19.fa.gz -O hg19.dict
  177.  
  178. java -jar gatk-4.1.2.0/gatk-package-4.1.2.0-local.jar BaseRecalibrator -R hg19.fa -I $prefix-4c-GATKgroup.bam -known-sites All_20180423.vcf.gz -O $prefix-4c-GATKrecal.table
  179.  
  180. java -jar gatk-4.1.2.0/gatk-package-4.1.2.0-local.jar ApplyBQSR -R hg19.fa -I $prefix-4c-GATKrecal.table -I $prefix-4c-GATKgroup.bam -O $prefix-4c-GATKrecal.bam
  181.  
  182. java -jar gatk-4.1.2.0/gatk-package-4.1.2.0-local.jar HaplotypeCaller -R hg19.fa -I $prefix-4c-GATKrecal.bam --genotyping-mode DISCOVERY -O $prefix-4c-GATKcns.vcf
  183.  
  184. echo ""
  185. echo "=== Iniciando filtrado ==="
  186. echo ""
  187.  
  188. java -jar gatk-4.1.2.0/gatk-package-4.1.2.0-local.jar SelectVariants -R hg19.fa -V $prefix-4c-GATKcns.vcf -O $prefix-4c-GATKsnp.vcf -select-type SNP
  189.  
  190. java -jar gatk-4.1.2.0/gatk-package-4.1.2.0-local.jar VariantFiltration --filter-name "NO PASS" -filter "QD<2.0 || FS>60.0 || MQ<40.0|| MQRankSum<-12.5 ||ReadPosRankSum<-8.0" -R hg19.fa -V $prefix-4c-GATKsnp.vcf -O $prefix-4c-GATKsnp_filt.vcf
  191.  
  192. java -jar gatk-4.1.2.0/gatk-package-4.1.2.0-local.jar SelectVariants -R hg19.fa -V $prefix-4c-GATKcns.vcf -O $prefix-4c-GATKindel.vcf -select-type INDEL
  193.  
  194. java -jar gatk-4.1.2.0/gatk-package-4.1.2.0-local.jar VariantFiltration --filter-name "NO PASS" -filter "QD<2.0 || FS>200.0 || ReadPosRankSum<-20.0" -R hg19.fa -V $prefix-4c-GATKindel.vcf -O $prefix-4c-GATKindel_filt.vcf
  195.  
  196. echo ""
  197. echo "=== Combinando ==="
  198. echo ""
  199.  
  200.  
  201. gatk MergeVcfs -I $prefix-4c-GATKsnp_filt.vcf -I $prefix-4c-GATKindel_filt.vcf -R hg19.fa -O $prefix-4c-GATKmerge.vcf
  202.  
  203. gatk SelectVariants -R hg19.fa -V $prefix-4c-GATKmerge.vcf -O $prefix-4c-GATKmerge_pass.vcf -select 'vc.isNotFiltered()'
  204.  
  205. echo ""
  206. echo "=== Uniendo VCFs ==="
  207. echo ""
  208.  
  209. if test -f "GenomeAnalysisTK.jar"; then
  210.  
  211.  
  212. java -jar GenomeAnalysisTK.jar -T CombineVariants -R hg19.fa -minN 2 -V: $prefix-4a-BCF_filt.vcf -V: $prefix-4b-VSCANcns_filt.vcf -V: $prefix-4c-GATKmerge_pass.vcf -o $prefix-5_GATKunion.vcf
  213.  
  214. echo "=== VCFs unidos, saliendo... ==="
  215.  
  216. else
  217. echo ""
  218. echo "=== No se ha encontrado el archivo GenomeAnalysisTK.jar. Por favor descarguelo y ponga el jar en la carpeta ${HOME} ==="
  219. echo ""
  220.  
  221. fi
  222.  
  223.  
  224.  
  225. else
  226. echo "=== El archivo de entrada indicado no existe. Debe estar en la misma ruta del script ==="
  227. fi
  228. fi
Advertisement
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!