Perl & Genomics

Parserizzare un file VCF, fornito dall’utente, e stampare solo le righe del file che soddisfino questi criteri:

• Che abbiano un codice “rs” (ovvero sono annotate in dbSNP)
• La qualità sia superiore ad una soglia definita dall’utente
• Che siano SNP (Single Nucleotide Polimorphism)
• La percentuale di righe stampate (ovvero che soddifano i criteri suddetti) sul totale di righe presenti

Questo post contiene informazioni specifiche per il progetto hw14. La lista dei progetti è presente in questa pagina.

Formato VCF

Il formato VCF serve a conservare varianti (SNP, INDEL…) trovate in un progetto di risequenziamento. E’ un formato tabulare che ha delle righe di intestazione precedute da “#”. Ecco un esempio

##fileformat=VCFv4.0
##fileDate=20110705
##reference=1000GenomesPilot-NCBI37
##phasing=partial
##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples With Data">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##INFO=<ID=AF,Number=.,Type=Float,Description="Allele Frequency">
##INFO=<ID=AA,Number=1,Type=String,Description="Ancestral Allele">
##INFO=<id=db,number=0,type=flag,description="dbSNP membership, build 129">
##INFO=<id=H2,Number=0,Type=Flag,Description="HapMap2 membership">
##FILTER=<ID=q10,Description="Quality below 10">
##FILTER=<ID=s50,Description="Less than 50% of samples have data">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
##FORMAT=<ID=HQ,Number=2,Type=Integer,Description="Haplotype Quality">
#CHROM POS    ID        REF  ALT     QUAL FILTER INFO                              FORMAT      Sample1        Sample2        Sample3
2      4370   rs6057    G    A       29   .      NS=2;DP=13;AF=0.5;DB;H2           GT:GQ:DP:HQ 0|0:48:1:52,51 1|0:48:8:51,51 1/1:43:5:.,.
2      7330   .         T    A       3    q10    NS=5;DP=12;AF=0.017               GT:GQ:DP:HQ 0|0:46:3:58,50 0|1:3:5:65,3   0/0:41:3
2      110696 rs6055    A    G,T     67   PASS   NS=2;DP=10;AF=0.333,0.667;AA=T;DB GT:GQ:DP:HQ 1|2:21:6:23,27 2|1:2:0:18,2   2/2:35:4
2      130237 .         T    .       47   .      NS=2;DP=16;AA=T                   GT:GQ:DP:HQ 0|0:54:7:56,60 0|0:48:4:56,51 0/0:61:2
2      134567 microsat1 GTCT G,GTACT 50   PASS   NS=2;DP=9;AA=G                    GT:GQ:DP    0/1:35:4       0/2:17:2       1/1:40:3

Come si nota il file è tabulare (separato da TAB), e le prime colonne sono di facile interpretazione:

1. Chrom (cromosoma in cui compare la variante)
2. Pos (nucleotide in cui compare la variante)
3. ID (la variante è in database? se si si mette il suo ID. dbSNP ha id che iniziano per rs…)
4. Ref (che base o che basi ci sono nel reference in quella posizione?)
5. Alt (che alternative ho trovato  a quelle basi? Separate con virgola)
6. Qual (qualità)

Prendiamo l’ultima riga. Nel reference ho GTCT, nel mio progetto ho trovato due varianti alternative: G (ovvero una delezione di TCT) e GTACT (ovvero l’inserzione di una A dopo GT).

Note sul programma

Il programma deve stampare l’intestazione così com’è presente nel file originale (per le vostre prove potete usare quello riportato qui sopra, tratto da Wikipedia).

Poi valutare se le condizioni espresse nella consegna sono soddisfatte. Se si stampare tutta la riga, altrimenti saltarla.