VCF – Variant Calling Format 基因變異儲存格式

生物資訊領域是近幾年來很熱門的一個領域,本篇紀錄有關 VCF – Variant Calling Format 基因變異儲存格式檔案的處理過程,有關於 VCF 的簡單介紹可以參考連結,也可以參考  Wikipedia 裡面關於 VCF 的介紹,一個 VCF 檔案大致上長成以下的樣子:

##fileformat=VCFv4.3
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=file:///seq/references/1000GenomesPilot-NCBI36.fasta
##contig=<ID=20,length=62435964,assembly=B36,md5=f126cdf8a6e0c7f379d618ff66beb2da,species="Homo sapiens",taxonomy=x>
##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=A,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=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##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       NA00001         NA00002          NA00003
20     14370    rs6054257  G     A      29    PASS    NS=3;DP=14;AF=0.5;DB;H2           GT:GQ:DP:HQ  0|0:48:1:51,51  1|0:48:8:51,51   1/1:43:5:.,.
20     17330    .          T     A      3     q10     NS=3;DP=11;AF=0.017               GT:GQ:DP:HQ  0|0:49:3:58,50  0|1:3:5:65,3     0/0:41:3
20     1110696  rs6040355  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
20     1230237  .          T     .      47    PASS    NS=3;DP=13;AA=T                   GT:GQ:DP:HQ  0|0:54:7:56,60  0|0:48:4:51,51   0/0:61:2
20     1234567  microsat1  GTC   G,GTCT 50    PASS    NS=3;DP=9;AA=G                    GT:GQ:DP     0/1:35:4        0/2:17:2         1/1:40:3

REF:基因公版 Reference 上面的核苷酸序列。
ALT:sample 上面呈現的核苷酸,可能會有一個以上,或是一整個序列的變異。
Common Format 裡面通常包含以下的值(參考):

  • GT (Genotype) 這邊需要注意由於人類是兩倍體基因 (雙套染色體)
    • 0/0: homozygous reference.
    • 0/1: heterozygous, one from reference, one from alternate.
    • 1/1: homozygous alternate.
    • NA: 沒有訊息
  • AD (Allele Depth), DP (Filtered Depth) 是兩個重要的指標標示這個點位的深度,也就是在 Read Mapping 的時候實際上 Map 到的次數,有一些 Read Mapping 會有 Filtered 的設定造成 AD 與 DP 的差異!
    • AD: [4,0]
    • DP: 4
  • PL (Phred-scaled likelihoods of the possible genotypes)
    假設是對 ALT 只有一個值的情況之下,PL 會有三個值分別為 0/0, 0/1, 1/1 如果 ALT 有兩個讀值,則 PL 會有不同的呈現方式!
  • GQ (Quality of the assigned genotype)
    類似 PL 告知 GT 的可信賴度!

利用以下的資料作為例子,在 873762 的點位上,公版是 T,樣本觀察到變異為 G,GT 顯示等位基因呈現 heterogyzous 也就是說一個是 T 一個是 G,他們各自的 AD 為 173 和 141,合計 DP 為 282,應該是 173+141 之後再篩選掉一些,GQ 為 99 代表非常確定就是 0/1 的基因型,0/0, 0/1, 1/1 各自的正規化 likelihoods 值?

CHROM  POS    ID         REF ALT QUAL FILTER INFO FORMAT         NA12878 
1      873762 .          T   G   [CLIPPED]        GT:AD:DP:GQ:PL 0/1:173,141:282:99:255,0,255 
1      877664 rs3828047  A   G   [CLIPPED]        GT:AD:DP:GQ:PL 1/1:0,105:94:99:255,255,0 
1      899282 rs28548431 C   T   [CLIPPED]        GT:AD:DP:GQ:PL 0/1:1,3:4:26:103,0,26

如果想要修改 VCF 內部的內容的話應該要怎麼做呢?

方法一 (解壓縮修改後再重新壓縮):
gzip -d example.vcf.gz

此時會得到一個 example.vcf 檔案,可以透過 vim 修改,不過要特別注意要 Header 裡面的 Column Name 之間需要用 Tab 隔開在修改的時候要注意,接下來需要使用以下的指令進行壓縮(參考連結):

bgzip example.vcf

備註:這邊不能夠使用 gzip 或是 gunzip 等等的壓縮方式,否則會在讀取的時候會出現以下的錯誤 (java.util.zip.ZipException: File does not conform to block gzip format.):

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/local/lib/python3.7/dist-packages/hail/matrixtable.py", line 2428, in count
    return Env.backend().execute(ir)
  File "/usr/local/lib/python3.7/dist-packages/hail/backend/backend.py", line 208, in execute
    result = json.loads(self._jhc.backend().executeJSON(jir))
  File "/home/spark-current/python/lib/py4j-0.10.7-src.zip/py4j/java_gateway.py", line 1257, in __call__
  File "/usr/local/lib/python3.7/dist-packages/hail/utils/java.py", line 211, in deco
    'Error summary: %s' % (deepest, full, hail.__version__, deepest)) from None
hail.utils.java.FatalError: ZipException: File does not conform to block gzip format.

Java stack trace:
java.lang.RuntimeException: error while applying lowering 'InterpretNonCompilable'
	at is.hail.expr.ir.lowering.LoweringPipeline$anonfun$apply$1.apply(LoweringPipeline.scala:26)
	at is.hail.expr.ir.lowering.LoweringPipeline$anonfun$apply$1.apply(LoweringPipeline.scala:18)
	at scala.collection.IndexedSeqOptimized$class.foreach(IndexedSeqOptimized.scala:33)
	at scala.collection.mutable.WrappedArray.foreach(WrappedArray.scala:35)
	at is.hail.expr.ir.lowering.LoweringPipeline.apply(LoweringPipeline.scala:18)
	at is.hail.expr.ir.CompileAndEvaluate$._apply(CompileAndEvaluate.scala:28)
	at is.hail.backend.spark.SparkBackend.is$hail$backend$spark$SparkBackend$_execute(SparkBackend.scala:297)
	at is.hail.backend.spark.SparkBackend$anonfun$execute$1.apply(SparkBackend.scala:284)
	at is.hail.backend.spark.SparkBackend$anonfun$execute$1.apply(SparkBackend.scala:283)
	at is.hail.expr.ir.ExecuteContext$anonfun$scoped$1.apply(ExecuteContext.scala:15)
	at is.hail.expr.ir.ExecuteContext$anonfun$scoped$1.apply(ExecuteContext.scala:13)
	at is.hail.utils.package$.using(package.scala:604)
	at is.hail.annotations.Region$.scoped(Region.scala:18)
	at is.hail.expr.ir.ExecuteContext$.scoped(ExecuteContext.scala:13)
	at is.hail.expr.ir.ExecuteContext$.scoped(ExecuteContext.scala:10)
	at is.hail.backend.spark.SparkBackend.execute(SparkBackend.scala:283)
	at is.hail.backend.spark.SparkBackend.executeJSON(SparkBackend.scala:303)
	at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
	at sun.reflect.NativeMethodAccessorImpl.invoke(NativeMethodAccessorImpl.java:62)
	at sun.reflect.DelegatingMethodAccessorImpl.invoke(DelegatingMethodAccessorImpl.java:43)
	at java.lang.reflect.Method.invoke(Method.java:498)
	at py4j.reflection.MethodInvoker.invoke(MethodInvoker.java:244)
	at py4j.reflection.ReflectionEngine.invoke(ReflectionEngine.java:357)
	at py4j.Gateway.invoke(Gateway.java:282)
	at py4j.commands.AbstractCommand.invokeMethod(AbstractCommand.java:132)
	at py4j.commands.CallCommand.execute(CallCommand.java:79)
	at py4j.GatewayConnection.run(GatewayConnection.java:238)
	at java.lang.Thread.run(Thread.java:748)

org.apache.spark.SparkException: Job aborted due to stage failure: Task 1 in stage 35.0 failed 1 times, most recent failure: Lost task 1.0 in stage 35.0 (TID 1433, localhost, executor driver): java.util.zip.ZipException: File does not conform to block gzip format.
	at is.hail.io.compress.BGzipInputStream$BGzipHeader.<init>(BGzipInputStream.java:36)
	at is.hail.io.compress.BGzipInputStream.decompressNextBlock(BGzipInputStream.java:150)
	at is.hail.io.compress.BGzipInputStream.<init>(BGzipInputStream.java:111)
	at is.hail.io.compress.BGzipCodec.createInputStream(BGzipCodec.java:46)
	at org.apache.hadoop.mapred.LineRecordReader.<init>(LineRecordReader.java:113)
	at org.apache.hadoop.mapred.TextInputFormat.getRecordReader(TextInputFormat.java:67)
	at org.apache.spark.rdd.HadoopRDD$anon$1.liftedTree1$1(HadoopRDD.scala:267)
	at org.apache.spark.rdd.HadoopRDD$anon$1.<init>(HadoopRDD.scala:266)
	at org.apache.spark.rdd.HadoopRDD.compute(HadoopRDD.scala:224)
	at org.apache.spark.rdd.HadoopRDD.compute(HadoopRDD.scala:95)
	at org.apache.spark.rdd.RDD.computeOrReadCheckpoint(RDD.scala:324)
	at org.apache.spark.rdd.RDD.iterator(RDD.scala:288)
	at org.apache.spark.rdd.MapPartitionsRDD.compute(MapPartitionsRDD.scala:52)
	at org.apache.spark.rdd.RDD.computeOrReadCheckpoint(RDD.scala:324)
	at org.apache.spark.rdd.RDD.iterator(RDD.scala:288)
	at org.apache.spark.rdd.MapPartitionsRDD.compute(MapPartitionsRDD.scala:52)
	at org.apache.spark.rdd.RDD.computeOrReadCheckpoint(RDD.scala:324)
	at org.apache.spark.rdd.RDD.iterator(RDD.scala:288)
	at org.apache.spark.rdd.MapPartitionsRDD.compute(MapPartitionsRDD.scala:52)
	at org.apache.spark.rdd.RDD.computeOrReadCheckpoint(RDD.scala:324)
	at org.apache.spark.rdd.RDD.iterator(RDD.scala:288)
	at org.apache.spark.rdd.MapPartitionsRDD.compute(MapPartitionsRDD.scala:52)
	at org.apache.spark.rdd.RDD.computeOrReadCheckpoint(RDD.scala:324)
	at org.apache.spark.rdd.RDD.iterator(RDD.scala:288)
	at org.apache.spark.rdd.MapPartitionsRDD.compute(MapPartitionsRDD.scala:52)
	at org.apache.spark.rdd.RDD.computeOrReadCheckpoint(RDD.scala:324)
	at org.apache.spark.rdd.RDD.iterator(RDD.scala:288)
	at org.apache.spark.scheduler.ResultTask.runTask(ResultTask.scala:90)
	at org.apache.spark.scheduler.Task.run(Task.scala:121)
	at org.apache.spark.executor.Executor$TaskRunner$anonfun$10.apply(Executor.scala:408)
	at org.apache.spark.util.Utils$.tryWithSafeFinally(Utils.scala:1360)
	at org.apache.spark.executor.Executor$TaskRunner.run(Executor.scala:414)
	at java.util.concurrent.ThreadPoolExecutor.runWorker(ThreadPoolExecutor.java:1149)
	at java.util.concurrent.ThreadPoolExecutor$Worker.run(ThreadPoolExecutor.java:624)
	at java.lang.Thread.run(Thread.java:748)

Driver stacktrace:
	at org.apache.spark.scheduler.DAGScheduler.org$apache$spark$scheduler$DAGScheduler$failJobAndIndependentStages(DAGScheduler.scala:1889)
	at org.apache.spark.scheduler.DAGScheduler$anonfun$abortStage$1.apply(DAGScheduler.scala:1877)
	at org.apache.spark.scheduler.DAGScheduler$anonfun$abortStage$1.apply(DAGScheduler.scala:1876)
	at scala.collection.mutable.ResizableArray$class.foreach(ResizableArray.scala:59)
	at scala.collection.mutable.ArrayBuffer.foreach(ArrayBuffer.scala:48)
	at org.apache.spark.scheduler.DAGScheduler.abortStage(DAGScheduler.scala:1876)
	at org.apache.spark.scheduler.DAGScheduler$anonfun$handleTaskSetFailed$1.apply(DAGScheduler.scala:926)
	at org.apache.spark.scheduler.DAGScheduler$anonfun$handleTaskSetFailed$1.apply(DAGScheduler.scala:926)
	at scala.Option.foreach(Option.scala:257)
	at org.apache.spark.scheduler.DAGScheduler.handleTaskSetFailed(DAGScheduler.scala:926)
	at org.apache.spark.scheduler.DAGSchedulerEventProcessLoop.doOnReceive(DAGScheduler.scala:2110)
	at org.apache.spark.scheduler.DAGSchedulerEventProcessLoop.onReceive(DAGScheduler.scala:2059)
	at org.apache.spark.scheduler.DAGSchedulerEventProcessLoop.onReceive(DAGScheduler.scala:2048)
	at org.apache.spark.util.EventLoop$anon$1.run(EventLoop.scala:49)
	at org.apache.spark.scheduler.DAGScheduler.runJob(DAGScheduler.scala:737)
	at org.apache.spark.SparkContext.runJob(SparkContext.scala:2061)
	at org.apache.spark.SparkContext.runJob(SparkContext.scala:2082)
	at org.apache.spark.SparkContext.runJob(SparkContext.scala:2114)
	at is.hail.sparkextras.ContextRDD.crunJobWithIndex(ContextRDD.scala:228)
	at is.hail.rvd.RVD$.getKeyInfo(RVD.scala:1226)
	at is.hail.rvd.RVD$.makeCoercer(RVD.scala:1301)
	at is.hail.io.vcf.MatrixVCFReader.coercer(LoadVCF.scala:1565)
	at is.hail.io.vcf.MatrixVCFReader.apply(LoadVCF.scala:1592)
	at is.hail.expr.ir.TableRead.execute(TableIR.scala:312)
	at is.hail.expr.ir.Interpret$anonfun$run$1.apply$mcJ$sp(Interpret.scala:598)
	at is.hail.expr.ir.Interpret$anonfun$run$1.apply(Interpret.scala:598)
	at is.hail.expr.ir.Interpret$anonfun$run$1.apply(Interpret.scala:598)
	at scala.Option.getOrElse(Option.scala:121)
	at is.hail.expr.ir.Interpret$.run(Interpret.scala:598)
	at is.hail.expr.ir.Interpret$.alreadyLowered(Interpret.scala:53)
	at is.hail.expr.ir.InterpretNonCompilable$.interpretAndCoerce$1(InterpretNonCompilable.scala:16)
	at is.hail.expr.ir.InterpretNonCompilable$.is$hail$expr$ir$InterpretNonCompilable$rewrite$1(InterpretNonCompilable.scala:53)
	at is.hail.expr.ir.InterpretNonCompilable$anonfun$1.apply(InterpretNonCompilable.scala:25)
	at is.hail.expr.ir.InterpretNonCompilable$anonfun$1.apply(InterpretNonCompilable.scala:25)
	at scala.collection.TraversableLike$anonfun$map$1.apply(TraversableLike.scala:234)
	at scala.collection.TraversableLike$anonfun$map$1.apply(TraversableLike.scala:234)
	at scala.collection.mutable.ResizableArray$class.foreach(ResizableArray.scala:59)
	at scala.collection.mutable.ArrayBuffer.foreach(ArrayBuffer.scala:48)
	at scala.collection.TraversableLike$class.map(TraversableLike.scala:234)
	at scala.collection.AbstractTraversable.map(Traversable.scala:104)
	at is.hail.expr.ir.InterpretNonCompilable$.rewriteChildren$1(InterpretNonCompilable.scala:25)
	at is.hail.expr.ir.InterpretNonCompilable$.is$hail$expr$ir$InterpretNonCompilable$rewrite$1(InterpretNonCompilable.scala:54)
	at is.hail.expr.ir.InterpretNonCompilable$.apply(InterpretNonCompilable.scala:58)
	at is.hail.expr.ir.lowering.InterpretNonCompilablePass$.transform(LoweringPass.scala:50)
	at is.hail.expr.ir.lowering.LoweringPass$anonfun$apply$3$anonfun$1.apply(LoweringPass.scala:15)
	at is.hail.expr.ir.lowering.LoweringPass$anonfun$apply$3$anonfun$1.apply(LoweringPass.scala:15)
	at is.hail.utils.ExecutionTimer.time(ExecutionTimer.scala:69)
	at is.hail.expr.ir.lowering.LoweringPass$anonfun$apply$3.apply(LoweringPass.scala:15)
	at is.hail.expr.ir.lowering.LoweringPass$anonfun$apply$3.apply(LoweringPass.scala:13)
	at is.hail.utils.ExecutionTimer.time(ExecutionTimer.scala:69)
	at is.hail.expr.ir.lowering.LoweringPass$class.apply(LoweringPass.scala:13)
	at is.hail.expr.ir.lowering.InterpretNonCompilablePass$.apply(LoweringPass.scala:45)
	at is.hail.expr.ir.lowering.LoweringPipeline$anonfun$apply$1.apply(LoweringPipeline.scala:20)
	at is.hail.expr.ir.lowering.LoweringPipeline$anonfun$apply$1.apply(LoweringPipeline.scala:18)
	at scala.collection.IndexedSeqOptimized$class.foreach(IndexedSeqOptimized.scala:33)
	at scala.collection.mutable.WrappedArray.foreach(WrappedArray.scala:35)
	at is.hail.expr.ir.lowering.LoweringPipeline.apply(LoweringPipeline.scala:18)
	at is.hail.expr.ir.CompileAndEvaluate$._apply(CompileAndEvaluate.scala:28)
	at is.hail.backend.spark.SparkBackend.is$hail$backend$spark$SparkBackend$_execute(SparkBackend.scala:297)
	at is.hail.backend.spark.SparkBackend$anonfun$execute$1.apply(SparkBackend.scala:284)
	at is.hail.backend.spark.SparkBackend$anonfun$execute$1.apply(SparkBackend.scala:283)
	at is.hail.expr.ir.ExecuteContext$anonfun$scoped$1.apply(ExecuteContext.scala:15)
	at is.hail.expr.ir.ExecuteContext$anonfun$scoped$1.apply(ExecuteContext.scala:13)
	at is.hail.utils.package$.using(package.scala:604)
	at is.hail.annotations.Region$.scoped(Region.scala:18)
	at is.hail.expr.ir.ExecuteContext$.scoped(ExecuteContext.scala:13)
	at is.hail.expr.ir.ExecuteContext$.scoped(ExecuteContext.scala:10)
	at is.hail.backend.spark.SparkBackend.execute(SparkBackend.scala:283)
	at is.hail.backend.spark.SparkBackend.executeJSON(SparkBackend.scala:303)
	at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
	at sun.reflect.NativeMethodAccessorImpl.invoke(NativeMethodAccessorImpl.java:62)
	at sun.reflect.DelegatingMethodAccessorImpl.invoke(DelegatingMethodAccessorImpl.java:43)
	at java.lang.reflect.Method.invoke(Method.java:498)
	at py4j.reflection.MethodInvoker.invoke(MethodInvoker.java:244)
	at py4j.reflection.ReflectionEngine.invoke(ReflectionEngine.java:357)
	at py4j.Gateway.invoke(Gateway.java:282)
	at py4j.commands.AbstractCommand.invokeMethod(AbstractCommand.java:132)
	at py4j.commands.CallCommand.execute(CallCommand.java:79)
	at py4j.GatewayConnection.run(GatewayConnection.java:238)
	at java.lang.Thread.run(Thread.java:748)

java.util.zip.ZipException: File does not conform to block gzip format.
	at is.hail.io.compress.BGzipInputStream$BGzipHeader.<init>(BGzipInputStream.java:36)
	at is.hail.io.compress.BGzipInputStream.decompressNextBlock(BGzipInputStream.java:150)
	at is.hail.io.compress.BGzipInputStream.<init>(BGzipInputStream.java:111)
	at is.hail.io.compress.BGzipCodec.createInputStream(BGzipCodec.java:46)
	at org.apache.hadoop.mapred.LineRecordReader.<init>(LineRecordReader.java:113)
	at org.apache.hadoop.mapred.TextInputFormat.getRecordReader(TextInputFormat.java:67)
	at org.apache.spark.rdd.HadoopRDD$anon$1.liftedTree1$1(HadoopRDD.scala:267)
	at org.apache.spark.rdd.HadoopRDD$anon$1.<init>(HadoopRDD.scala:266)
	at org.apache.spark.rdd.HadoopRDD.compute(HadoopRDD.scala:224)
	at org.apache.spark.rdd.HadoopRDD.compute(HadoopRDD.scala:95)
	at org.apache.spark.rdd.RDD.computeOrReadCheckpoint(RDD.scala:324)
	at org.apache.spark.rdd.RDD.iterator(RDD.scala:288)
	at org.apache.spark.rdd.MapPartitionsRDD.compute(MapPartitionsRDD.scala:52)
	at org.apache.spark.rdd.RDD.computeOrReadCheckpoint(RDD.scala:324)
	at org.apache.spark.rdd.RDD.iterator(RDD.scala:288)
	at org.apache.spark.rdd.MapPartitionsRDD.compute(MapPartitionsRDD.scala:52)
	at org.apache.spark.rdd.RDD.computeOrReadCheckpoint(RDD.scala:324)
	at org.apache.spark.rdd.RDD.iterator(RDD.scala:288)
	at org.apache.spark.rdd.MapPartitionsRDD.compute(MapPartitionsRDD.scala:52)
	at org.apache.spark.rdd.RDD.computeOrReadCheckpoint(RDD.scala:324)
	at org.apache.spark.rdd.RDD.iterator(RDD.scala:288)
	at org.apache.spark.rdd.MapPartitionsRDD.compute(MapPartitionsRDD.scala:52)
	at org.apache.spark.rdd.RDD.computeOrReadCheckpoint(RDD.scala:324)
	at org.apache.spark.rdd.RDD.iterator(RDD.scala:288)
	at org.apache.spark.rdd.MapPartitionsRDD.compute(MapPartitionsRDD.scala:52)
	at org.apache.spark.rdd.RDD.computeOrReadCheckpoint(RDD.scala:324)
	at org.apache.spark.rdd.RDD.iterator(RDD.scala:288)
	at org.apache.spark.scheduler.ResultTask.runTask(ResultTask.scala:90)
	at org.apache.spark.scheduler.Task.run(Task.scala:121)
	at org.apache.spark.executor.Executor$TaskRunner$anonfun$10.apply(Executor.scala:408)
	at org.apache.spark.util.Utils$.tryWithSafeFinally(Utils.scala:1360)
	at org.apache.spark.executor.Executor$TaskRunner.run(Executor.scala:414)
	at java.util.concurrent.ThreadPoolExecutor.runWorker(ThreadPoolExecutor.java:1149)
	at java.util.concurrent.ThreadPoolExecutor$Worker.run(ThreadPoolExecutor.java:624)
	at java.lang.Thread.run(Thread.java:748)
方法二 (一氣呵成處理完):

概念與方法一類似,但是把所有的步驟用一個指令貫穿完成,例如以下的指令:

zcat $f | sed 's:^chr::g' | sed 's:^M:MT:g' | bgzip > $2$fileName

如果想要批次處理一個目標資料夾裡面的所有 .vcf.gz 檔的話可以參考以下的做法。

FileFolder=$1
DestFolder=$2

for f in $FileFolder;
do
    echo "Start to proceed: "$f
    fileName=$(basename $f suffix)
    zcat $f | sed 's:^chr::g' | sed 's:^M:MT:g' | bgzip > $2$fileName
done
gVCF – Genomic VCF (參考連結)

gVCF 與 VCF 最主要的差異在於 VCF 只記錄有變異的基因點位(variant site record),gVCF 則是會另外記錄沒有變異的點位(non-var block record),這樣做的好處是在分析 cohort 也就是一群樣本的時候比較好做 join 的運算!一般來說 gVCF 是由 GATK version 3.0 中的 HaplotypeCaller 生產出來的!