1. fasta和fastq
1.1. fasta：序列
- 以 > 開頭
- gi|gi號|來源標識|序列標識（接收號/名稱等）
- 序列描述
- 序列中允許空格婉支，換行留搔，空行，直到下一個 > 表示序列結(jié)束
1.2. fastq：序列+質(zhì)量
- 第一行以 @ 開頭，唯一的序列ID標識符以现，可選的序列描述內(nèi)容，用空格分開
- 第二行：序列字符（核酸/氨基酸）
- 第三行：“+”開頭约啊，后面空著或者加上第一行@之后的內(nèi)容
- 第四行：堿基質(zhì)量字符邑遏，每個字符對應(yīng)第二行相應(yīng)位置堿基或氨基酸的質(zhì)量，可以按一定規(guī)則轉(zhuǎn)換為堿基質(zhì)量得分恰矩，進而反映該堿基的錯誤率记盒，這一行字符數(shù)與第二行字符數(shù)必須相等
1.3. 處理
-
fastqc test.fq即可得到html格式的報告，可以可視化地查看堿基是否平衡外傅，測序質(zhì)量如何等等
1.4. 作業(yè)
echo $(wc -l reads_1.fq|awk '{print $1}')/4|bc-
awk '/^@/{print $0}' reads_1.fq或awk '{if(NR%4==1)print}' read_1.fq cat reads_1.fq |paste - - - - |cut -f2cat reads_1.fq |paste - - - - |cut -f3cat reads_1.fq |paste - - - - |cut -f4cat reads_1.fq |paste - - - - |cut -f2|grep 'N'|wc-
cat reads_1.fq |paste - - - - |cut -f2|wc -m- 用上面這種方法會多出10000個回車符纪吮，所以要么手動減去10000俩檬，要么用下面的方法
awk '{if(NR%4==2)print length}' read_1.fq|paste -sd +|bccat reads_1.fq |paste - - - - |cut -f2|grep -o [ATCGN]|wc
cat reads_1.fq |paste - - - - |cut -f2|grep -o 'N'|wccat reads_1.fq |paste - - - - |cut -f4|grep -o '5'|wccat reads_1.fq |paste - - - - |cut -f4|grep -o '?'|wc-
cat reads_1.fq |paste - - - - |cut -f2|grep '^[ATCGatcg]'|wc- 理解錯了題目意思，正解：
awk '{if(NR%4==2)print length}' read_1.fq|cut -c1|sort|uniq -c
- 理解錯了題目意思，正解：
cat reads_1.fq |paste - - - - |cut -f1,2|tr '\t' '\n'|tr '@' '>'>reads_1.faecho $(wc -l reads_1.fa|awk '{print $1}')/2|bccat reads_1.fa |paste - - |cut -f2|grep -o [CG]|wccat reads_1.fa |tr -d 'N'cat reads_1.fa |grep -v 'N'-
cat reads_1.fa |paste - - |cut -f2|awk '{if(length>=65)print}'- 這樣會把第一行的序列信息去掉碾盟，用下面的方法更好
cat reads_1.fa |paste - - |awk '{if(length($2)>=65)print}'
-
awk '{if(NR%2==0){print substr($0,6,length($0)-10)}}' reads_1.fa這個做法的缺點同樣是丟失了序列信息 cat reads_1.fa |paste - - |awk '{if(length($2)<=125)print}'- 直接用fastqc
測序質(zhì)量.png
2. sam和bam
2.1. 標頭注釋部分
- 標頭以@開頭棚辽，用不同的tag表示不同的信息，主要有：
- @HQ冰肴，說明符合標準的版本屈藐、對比序列的排列順序
- @SQ，參考序列說明
- @PG熙尉，使用的比對程序說明
- LN联逻，參考序列的長度
2.2. 比對結(jié)果部分
- 每一行表示一個read的比對信息
- 每行包括11個必須字段和1個可選字段，字段之間用tag分割
2.3. 作業(yè)
cat tmp.sam |grep -v '@'|cut -f1|uniq|wccat tmp.sam |grep -v '@'|cut -f2|sort|uniq -c|sort -k1,1 -ncat tmp.sam |grep -v '^@' |awk '{if($6=="*")print$10}'|head -3
cat tmp.sam |grep -v '^@' |awk '{if($6=="*")print$1}'|sort|uniq -c|grep -w 1
cat tmp.sam |grep -v '^@' |awk '{if($6=="*")print$1}'|sort|uniq -c|grep -w 2
cat tmp.sam |grep -v '^@' |awk '{if($5>30)print}'cat tmp.sam |grep -v '^@' |awk '{if($6!="*")print$6}'|grep '[IDNSHPX]'cat tmp.sam |grep -v '^@' |awk '{if($9>1250||$9<-1250)print}'cat tmp.sam |grep -v '^@' |cut -f3 |sort|uniq -ccat tmp.sam |grep -v '^@' |awk '{if($10 ~/N/) print $0}'awk '{if($10 ~/N/) print}' tmp.sam |awk '{if($6 !~ /[IDNSHPX*]/) print}'grep '^@' tmp.sam |wcgrep -v '^@' tmp.sam |cut -f 12-1000|head -6grep -v '^@' tmp.sam |cut -f 12-1000|awk '{print NF}' |sort |uniq -cgrep '^@SQ' tmp.sam |awk '{print $3}'|cut -d: -f2awk '{if($6 ~/I/)print}' tmp.samawk '{if($6 ~/D/)print}' tmp.samawk '{if($4>5013 && $4<50130) print}' tmp.samcat tmp.sam |grep -v '^@' |awk '{if($4!=0)print}'|sort -nk4,4|less -SNgrep "102M3D11M" tmp.sam |awk '{print length($10)}'
2.4. 其它概念題
- 高通量測序數(shù)據(jù)分析中检痰，序列與參考序列的比對后得到的標準文件為什么文件包归？
A：sam文件
- sam文件是一種比對后的文件，能直接查看嗎攀细？
A：可以
- Sam/Bam文件分為兩部分箫踩，一部分以@開頭的是什么序列，另一部分是什么序列谭贪？
A：分別是標頭注釋部分和比對結(jié)果部分
- Sam文件除頭文件境钟，以什么符號分割文本的，比對部分信息一行是多少列俭识？你能用awk計算列數(shù)嗎慨削？
A：以tab分割，必須字段11列套媚，可選字段不一定缚态，可以用awk '{print NF}'
- Sam/Bam文件的@SQ開頭的行是什么？你能生成一個文本堤瘤，兩列玫芦，一列是參考基因組染色體/sca id，一列是長度嗎本辐？
A：參考染色體的信息桥帆，cat tmp.sam|grep -w ^@SQ|cut -f2-3 > chrom.txt
- Sam文件的比對位置是從1還是0開始的？
A：從1開始
- 常見CIGAR字符串各字母代表的意義
A：
- M：alignment match (can be a sequence match or mismatch)表示read可mapping到第三列的序列上慎皱，則read的堿基序列與第三列的序列堿基相同老虫，表示正常的mapping結(jié)果，M表示完全匹配茫多，但是無論reads與序列的正確匹配或是錯誤匹配該位置都顯示為MI：insertion to the reference表示read的堿基序列相對于第三列的RNAME序列祈匙，有堿基的插入
- D：deletion from the reference表示read的堿基序列相對于第三列的RNAME序列，有堿基的刪除
- N：skipped region from the reference表示可變剪接位置
- P：padding (silent deletion from padded reference)
- S：soft clipping (clipped sequences present in SEQ)
- H：hard clipping (clipped sequences NOT present in SEQ)clipped均表示一條read的序列被分開，之所以被分開夺欲，是因為read的一部分序列能匹配到第三列的RNAME序列上跪帝，而被分開的那部分不能匹配到RNAME序列上。
- “=”表示正確匹配到序列上
- “X”表示錯誤匹配到序列上
- 比對質(zhì)量的數(shù)字是哪一列些阅？越大比對質(zhì)量越好還是越小越好歉甚？
A：第五列，越大越好扑眉。
- Sam文件的flag是第幾列，數(shù)字代表什么意義赖钞？是怎么計算來的腰素？
A：第二列
- 1（1）該read是成對的paired reads中的一個
- 2（10）paired reads中每個都正確比對到參考序列上
- 4（100）該read沒比對到參考序列上
- 8（1000）與該read成對的matepair read沒有比對到參考序列上
- 16（10000）該read其反向互補序列能夠比對到參考序列
- 32（100000）與該read成對的matepair read其反向互補序列能夠比對到參考序列
- 64（1000000）在paired reads中，該read是與參考序列比對的第一條
- 128（10000000）在paired reads中雪营，該read是與參考序列比對的第二條
- 256（100000000）該read是次優(yōu)的比對結(jié)果
- 512（1000000000）該read沒有通過質(zhì)量控制
- 1024（10000000000）由于PCR或測序錯誤產(chǎn)生的重復(fù)reads
- 2048（100000000000）補充匹配的read
多個數(shù)字相加得到
- Sam文件怎么轉(zhuǎn)bam文件弓千？用什么指令？為什么要轉(zhuǎn)換献起？
A：samtools view -bS tmp.sam >tmp.bam 因為bam占用的磁盤空間比sam文本文件醒蠓谩；利用bam二進制文件的運算速度快
- Bam文件查看用什么指令谴餐？如果需要查看頭文件需要增加參數(shù)姻政？
A：samtools view tmp.bam，加上選項-H
- Bam文件為什么要排序岂嗓？排序后的bam和未排序的bam頭文件和chr position列有什么區(qū)別汁展？
A：排序后所有reads按照其在參考染色體上的位置排好
- Bam文件建索引的指令是什么？
A：samtools index <in.bam> [out.index]
- Bam文件可視化用什么工具厌殉？查看時需要建立索引嗎食绿？
A：igv 需要
- Bam文件統(tǒng)計reads比對情況用samtools的哪一個子命令？
A：samtools idxstats tmp.bam
- SE測序和PE測序的所比對后得到的sam文件的區(qū)別在哪里公罕？
A：qname第一列器紧，雙端有重復(fù)現(xiàn)象，單端沒有重復(fù)現(xiàn)象
- Bam能用gzip再壓縮嗎楼眷？
A：不能铲汪，雖然可以執(zhí)行tar -zcvf tmp.tar.gz tmp.bam命令，但隨后用ls -lh查看可以發(fā)現(xiàn)文件大小沒有變化
結(jié)果示意.png
- Sam文件通常由哪些比對軟件得到的摩桶？
A：bowtie2桥状，bwa
- Sam/Bam文件能轉(zhuǎn)成fastq文件嗎？
A：能
samtools view -bS aln.sam > aln.bam
samtools view -h -o aln.sam aln.bam
bam2fastq --aligned input.bam -o output.fq
- 有時候不能通過文件名的后綴來區(qū)別是sam還是bam硝清，你是怎么區(qū)別的辅斟？
A：能否直接cat查看
3. gff和gtf
- gff主要用來注釋基因組
- gtf主要用來注釋基因
3.1. gtf
- 每一列的含義
- seqid 參考序列ID
- source 注釋的來源 未知則用 . 代替。一般指明產(chǎn)生此文件的軟件或方法
- type 類型 此處名詞相對自由芦拿，建議使用符合SO管理的名稱士飒，如gene查邢，repeat_region，exon酵幕，CDS等
- start 開始位點 從1開始計數(shù)
- end 結(jié)束位點 對于一些可以量化的屬性扰藕，可以在此設(shè)置一個數(shù)值以表示程度的不同。如果為空芳撒，用 . 代替邓深。
- score 得分
- strand 正/負鏈
- phase 步進 對于編碼蛋白質(zhì)的CDS，本列指定下一個密碼子開始的位置笔刹，可以是0芥备、1或2，表示到達下一個密碼子要跳過的堿基個數(shù)
- attributes 屬性 格式為標簽=值（tag=value）不同屬性之間以分號相隔
3.2. gff
- 每一列的含義
- seqname 序列的名字 通常格式染色體ID或是contig ID
- source 注釋的來源 通常是預(yù)測軟件名或是公共數(shù)據(jù)庫
- start 開始位點 從1開始計數(shù)
- end 結(jié)束位點
- feature 基因結(jié)構(gòu) CDS舌菜，start_condon萌壳，stop_condon是一定要含有的類型
- score 得分 對該類型存在性和其坐標的可信度，非必須日月，可以用點 . 代替
- strand 正/負鏈 鏈的正向與負向袱瓮，分別用+/-表示
- frame 密碼子偏移 可以是0，1或2
- attributes 屬性
- gene_id value：表示轉(zhuǎn)錄本在基因組上的基因座的唯一的ID爱咬。
- gene_id與value值用空格分開尺借，如果值為空，則表示沒有對應(yīng)的基因精拟。
- transcript_id value：預(yù)測的轉(zhuǎn)錄本的唯一ID
- transcript_id與value值用空格分開褐望，空表示沒有轉(zhuǎn)錄本。
3.3. 關(guān)系
- gtf2和gff3的內(nèi)容很相似串前，區(qū)別只在其中的兩列
| gtf2 | gff3 | |
|---|---|---|
| feature/type | 必須注明 | 可以是任意名稱 |
| attributes | 名稱和值以空格隔開 | 名稱和值以=隔開 |
- 兩種文件格式之間的轉(zhuǎn)換：使用Cufflinks里的工具“gffread”
#gff2gtf
gffread my.gff3 -T -o my.gtf
#gtf2gff
gffread merged.gtf -o- >merged.gff3
4. vcf
- VCF格式瘫里，Variant Call Format 變異判讀文件格式
- 分為兩部分內(nèi)容：以“#”開頭的注釋部分；沒有“#”開頭的主體部分
- VCF文件主題部分的結(jié)構(gòu)
- CHROM：表示變異位點是在哪個contig里call出來的荡碾，如果是人類全基因組的話那就是chr1…chr22, chrX, Y, M
- POS： 變異位點相對于參考基因組所在的位置谨读，如果是indel，就是第一個堿基所在的位置
- ID： 如果call出來的SNP存在于dbsnp數(shù)據(jù)庫里坛吁，就會顯示相應(yīng)的dbsnp里的rs編號, 若沒有劳殖，則用 . 表示其為一個novel variant。
- REF： 在這個變異位點處拨脉，參考基因組中所對應(yīng)的堿基
- ALT： 在這個變異位點處哆姻，研究對象基因組中所對應(yīng)的堿基
- QUAL：Phred格式(Phred_scaled)的質(zhì)量值，表示在該位點存在variant的可能性玫膀；該值越高矛缨，則variant的可能性越大；計算方法：Phred值 = -10 * l/g (1-p) p為variant存在的概率; 通過計算公式可以看出值為10的表示錯誤概率為0.1，該位點為variant的概率為90%箕昭。
- FILTER：使用上一個QUAL值來進行過濾的話灵妨，是不夠的。GATK能使用其它的方法來進行過濾落竹，過濾結(jié)果中通過則該值為'PASS'泌霍；若variant不可靠，則該項不為'PASS'或'.'述召。
- INFO
- FORMAT 和 testxxx：這兩行合起來提供了'testxxx'這個sample的基因型的信息朱转。 testxxx 代表這該名稱的樣品，是由BAM文件中的@RG下的 SM 標簽決定的积暖。
- GT：樣品的基因型（genotype）肋拔。兩個數(shù)字中間用'/'分開，這兩個數(shù)字表示雙倍體的sample的基因型呀酸。0 表示樣品中有ref的allele；1 表示樣品中variant的allele琼梆； 2表示有第二個variant的allele性誉。因此：0/0 表示sample中該位點為純合的，和ref一致茎杂；0/1 表示sample中該位點為雜合的错览，有ref和variant兩個基因型；1/1 表示sample中該位點為純合的煌往，和variant一致倾哺。
- AD：對應(yīng)兩個以逗號隔開的值，這兩個值分別表示覆蓋到REF和ALT堿基的reads數(shù)刽脖，相當(dāng)于支持REF和支持ALT的測序深度羞海。
- DP：覆蓋到這個位點的總的reads數(shù)量，相當(dāng)于這個位點的深度（并不是多有的reads數(shù)量曲管，而是大概一定質(zhì)量值要求的reads數(shù)）却邓。
- PL：指定的三種基因型的質(zhì)量值(provieds the likelihoods of the given genotypes)。這三種指定的基因型為(0/0,0/1,1/1)院水，這三種基因型的概率總和為1腊徙。和之前不一致，該值越大檬某，表明為該種基因型的可能性越小撬腾。 Phred值 = -10 * lg § p為基因型存在的概率=10^（-Phred值/10）。
INFO - AC：表示該Allele的數(shù)目恢恼，Allele數(shù)目為1表示雙倍體的樣本在該位點只有1個等位基因發(fā)生了突變
- AF：表示Allele的頻率民傻，Allele頻率為0.5表示雙倍體的樣本在該位點只有50%的等位基因發(fā)生了突變
- AN：表示Allele的總數(shù)目,即：對于1個diploid sample而言：則基因型 0/1 表示sample為雜合子，Allele數(shù)為1(雙倍體的sample在該位點只有1個等位基因發(fā)生了突變)，Allele的頻率為0.5(雙倍體的 sample在該位點只有50%的等位基因發(fā)生了突變)饰潜，總的Allele為2初坠； 基因型 1/1 則表示sample為純合的，Allele數(shù)為2彭雾，Allele的頻率為1碟刺，總的Allele為2。
- DP：樣本在這個位置的reads覆蓋度薯酝，是一些reads被過濾掉后的覆蓋度（跟上面提到的DP類似）
- FS：使用Fisher's精確檢驗來檢測strand bias而得到的Fhred格式的p值半沽，值越小越好
- MQ：表示覆蓋序列質(zhì)量的均方值RMS Mapping Quality
4.1. 作業(yè)
cat tmp.vcf|grep -v "^#"|grep INDEL``cat tmp.vcf|grep -v "^#"|grep -v INDEL
(test) [root@instance-ptelbf6dreads]# cat tmp.vcf|grep -v "^#"|grep INDEL|cut -f8|cut -d';' -f4|cut -d'=' -f2|paste -sd +|bc
2831
(test) [root@instance-ptelbf6dreads]# cat tmp.vcf|grep -v "^#"|grep INDEL|wc -l
70
(test) [root@instance-ptelbf6dreads]# echo 2831/70|bc
40
#以上方法可能略笨，而且只能得到整數(shù)吴菠，下面附上別人寫的更好的方法
(test) [root@instance-ptelbf6dreads]# cat tmp.vcf|grep -v "^#"|grep INDEL|cut -f8|cut -d";" -f4|cut -d"=" -f2|awk '{ sum += $0; } END { print sum/NR }'
40.4429
(test) [root@instance-ptelbf6dreads]# cat tmp.vcf|grep -v "^#"|grep -v INDEL|cut -f8|cut -d";" -f1|cut -d"=" -f2|awk '{ sum += $0; } END { print sum/NR }'
37
cat tmp.vcf|grep -v "^#"|grep INDEL|awk '{if(length($4)<length($5))print}' # insertion
cat tmp.vcf|grep -v "^#"|grep INDEL|awk '{if(length($4)>length($5))print}' # deletion
-
cat tmp.vcf|grep -v "^#"|grep -v INDEL|awk '{print $4"-->"$5}'|sort|uniq -c|sort -k1,1nrSNP行的第4者填、第5行分別表示參考基因和研究對象的堿基類型 -
cat tmp.vcf|grep -v "^#"|cut -f10|cut -d: -f1|sort|uniq -c缺陷是沒有把整條取出來 - ``暫略
cat tmp.vcf|grep -v "^#"|awk '{if($6>30)print}'- ``暫略
- ``暫略
- IGV暫略
4.2. 其他思考題
- vcf的全稱是什么？是用來記錄什么信息的標準格式的文本做葵？
A：VCF文件全稱為Variant Call Format占哟，表示基因組的變異信息。
- 一般選用哪個指令查看vcf文件酿矢，為什么不用vim?
A：一般選用less -SN榨乎，因為vim沒有單行顯示，看得不清楚瘫筐。
- vcf文件以'##'開頭的是什么信息蜜暑？請認真查看這些信息。'#'開頭的是什么信息策肝？
A：以##開頭的是mate-information line肛捍，以##開頭，格式為key=value之众。 fileformat是必須的字段拙毫，表明VCF格式的版本，其他行主要用來描述INFO, FORMAT, FILTER等字段的具體含義棺禾；以#開頭的是heaher line恬偷，只有一行，行內(nèi)以\t分隔帘睦，是下方具體信息的表頭袍患。
- vcf文件除頭信息，每行有多少列竣付，請詳細敘述每行的含義诡延！請準確記憶。
A：至少9列古胆。1）CHROM: 染色體名字肆良；2）POS: 染色體的位置筛璧，起始位置為1；3）ID: 變異位點在數(shù)據(jù)庫中的ID惹恃，如果是dbsnp數(shù)據(jù)庫夭谤，推薦使用rs號，如果沒有ID巫糙，用點號表示缺失值朗儒；4）REF: 參考基因組上的堿基；5）ALT: 變異之后的堿基参淹；6）QUAL: 變異位點的質(zhì)量醉锄，質(zhì)量值越高，為真實的變異位點的概率越大浙值；7）FILTER: 過濾信息恳不，PASS代表通過了過濾；對于過濾失敗的位點开呐，會給出對應(yīng)的過濾失敗的原因烟勋，具體的含義可以查看mate-information line 中對FILTER字段的描述；8）INFO：額外的信息筐付，具體的含義可以查看mate-information line 中對INFO 字段的描述卵惦，看上去是一列，但其中的內(nèi)容可以無限擴增家妆；9）Format：表頭的##FORMAT就是對第九列的解釋，主要包括某一個特定位點基因型冕茅、測序深度的表述伤极。
- 理解format列和樣本列的對應(yīng)關(guān)系以及GT AD DP的含義。
A：GT：genotype姨伤；AD：Allele Depth哨坪，樣本中每一種allel的reads覆蓋度，在二倍體中是用逗號分隔的兩個數(shù)乍楚，前面對應(yīng)ref当编，后面對應(yīng)variant；DP：Raw read depth
- vcf文件第三列如果不是'.'徒溪，出現(xiàn)的rs號的id是什么忿偷？
A：變異位點名稱，對應(yīng)dbSNP數(shù)據(jù)庫中的ID
7.** vcf文件的ref臊泌，alt列和樣本列的0/1 1/1 或者1/2的聯(lián)系鲤桥？**
A：0代表樣本中ref的allel；1代表樣本variant的allel渠概；2表示有第二個variant的allel茶凳；0/0表示樣本中該位點純合嫂拴，與ref一致；0/1表示樣本中該位點雜合贮喧，有ref和variant兩個基因型筒狠；1/1表示樣本中位點純合，與variant一致
- vcf文件一般用什么軟件生成箱沦？請至少說出兩個軟件辩恼。請注意不同軟件生成的vcf格式的稍有不同的地方。
A：GATK和Samtools
- Vcf文件一般都比較大饱普，壓縮后的.gz文件用什么指令直接查看而不用解壓后查看运挫？
A：zless zcat
- 了解gvcf是什么格式，gvcf全稱是什么套耕？他與vcf有什么前后聯(lián)系谁帕？
A：GVCF代表Genomic VCF，GVCF是一種VCF冯袍，因此基本格式規(guī)范與常規(guī)VCF相同匈挖，但基因組VCF包含額外信息。術(shù)語GVCF有時僅用于描述包含基因組中每個位置（或感興趣區(qū)間）的記錄的VCF康愤，無論是否在該位點檢測到變體（例如由UnifiedGenotyper產(chǎn)生的VCF --output_mode EMIT_ALL_SITES）儡循。HaplotypeCaller 3.x生成的GVCF包含以非常特定的方式格式化的其他信息。
- 把alt列出現(xiàn)','的行提取出來
A：cat tmp.vcf |grep -v "^#"|awk '{if($5 ~",")print}'
- 請將chrid征冷、postion择膝、ref、alt检激、format肴捉、樣本列切割出來生成一個文本
A：cat tmp.vcf |grep -v "^#"|awk '{print $1,$2,$4,$5,$9}'> tmp.new.txt
- 將一個含snp，indel信息的vcf拆成一個只含snp叔收，一個只含indel信息的2個vcf文件齿穗。可借鑒軟件
A：
cat tmp.vcf |grep -v "^#"|awk '{if($8~"INDEL")print}'>indel.vcf
cat tmp.vcf |grep -v "^#"|awk '{if($8!~"INDEL")print}'>snp.vcf
- 用指令操作indel的vcf文件饺律，提取indel長度>4的變異行數(shù)窃页，存成一個文本。
A：cat indel.vcf |cut -f8|cut -d';' -f2無法提取完整信息
用vcftools過濾vcf文件复濒，如maf設(shè)置成0.05脖卖， depth設(shè)置成5-20，統(tǒng)計過濾前后的變異位點的總個數(shù)
利用vcftools提取每個樣本每一個位點的變異信息和深度信息巧颈，生成一個矩陣的文件胚嘲，至少含chrid，postion洛二，sample_DP馋劈，sample_GT四項信息
vcftools暫略
- 提取出變異位點上樣本有純和突變的行數(shù)
A：grep 1/1 snp.vcf |wc
- 統(tǒng)計一下1號染色體上的變異總個數(shù)攻锰。
A：示例文件的序列全部在一條染色體上。
- 提取一下BRCA1基因上發(fā)生變異的行數(shù)妓雾，如果是人類wes變異結(jié)果文件
暫略
- 統(tǒng)計vcf文件各樣本的缺失率娶吞，如果是多個樣本的群體call結(jié)果。
暫略
5. 小技巧
-
find . -name *.sam搜索sam格式文件
6. 本文引用源
最后械姻，向大家隆重推薦生信技能樹的一系列干貨妒蛇！
- 生信技能樹全球公益巡講：https://mp.weixin.qq.com/s/E9ykuIbc-2Ja9HOY0bn_6g
- B站公益74小時生信工程師教學(xué)視頻合輯https://mp.weixin.qq.com/s/IyFK7l_WBAiUgqQi8O7Hxw
- 招學(xué)徒https://mp.weixin.qq.com/s/KgbilzXnFjbKKunuw7NVfw
