首先使用bowtie2軟件自帶的測試數(shù)據(jù)生成sam/bam文件鸥咖，還有vcf文件
代碼如下：
mkdir -p ~/biosoft
cd ~/biosoft
wget https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh
bash Miniconda3-latest-Linux-x86_64.sh
source ~/.bashrc
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
conda config --set show_channel_urls yes
conda create -y -n test
conda activate test
conda install -y samtools bcftools
wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.3.4.3/bowtie2-2.3.4.3-linux-x86_64.zip
unzip bowtie2-2.3.4.3-linux-x86_64.zip
cd ~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads
../../bowtie2 -x ../index/lambda_virus -1 reads_1.fq -2 reads_2.fq | samtools sort -@ 5 -o tmp.bam
bcftools mpileup -f ../reference/lambda_virus.fa tmp.bam | bcftools call -vm > tmp.vcf
LINUX練習(xí)題
- 把突變記錄的vcf文件區(qū)分成 INDEL和SNP條目
grep 'INDEL' tmp.vcf | grep -v '#' # INDEL條目
grep -v -e 'INDEL' -v -e '#' tmp.vcf # SNP條目
- 統(tǒng)計INDEL和SNP條目的各自的平均測序深度
cat tmp.vcf | grep -v "^#" | grep 'INDEL' | cut -f 8 | cut -d ";" -f 4 | awk '{print substr($0, 4)}' | awk '{sum += $0} END {print sum/NR}' # INDEL條目
cat tmp.vcf | grep -v "^#" | grep -v 'INDEL' | cut -f 8 | cut -d ";" -f 1 | awk '{print substr($0, 4)}' | awk '{sum += $0} END {print sum/NR}' # SNP條目
- 把INDEL條目再區(qū)分成insertion和deletion情況
grep -v '^#' tmp.vcf | grep 'INDEL' | awk '{if ($4 > $5) print $0}' | less -SN # Deletion
grep -v '^#' tmp.vcf | grep 'INDEL' | awk '{if ($4 < $5) print $0}' | less -SN # Insertion
- 統(tǒng)計SNP條目的突變組合分布頻率
grep -v '^#' tmp.vcf | grep -v 'INDEL' | cut -f 4-5 | tr '\t' '-' | sort | uniq -c
# grep -v '^#' tmp.vcf | grep -v 'INDEL' | awk '{print $4"-"$5}'| sort | uniq -c
- 找到基因型不是 1/1 的條目燕鸽，個數(shù)
grep -v '^#' tmp.vcf | awk '{if ($10 !~ "1/1") print}'
grep -v '^#' tmp.vcf | awk '{if ($10 !~ "1/1") print}' | wc -l
- 篩選測序深度大于20的條目
egrep 'DP=2[1-9]|DP=[3-9][0-9]' tmp.vcf | less -SN
- 篩選變異位點質(zhì)量值大于30的條目
grep -v '^#' tmp.vcf | awk '{if ($6 > 30) print}' | less -SN
- 組合篩選變異位點質(zhì)量值大于30并且深度大于20的條目
grep -v '^#' tmp.vcf | awk '{if ($6 > 30) print}' | egrep 'DP=2[1-9]|DP=[3-9][0-9]' | less -SN
- 理解DP4=4,7,11,18 這樣的字段，就是 Number of high-quality ref-forward , ref-reverse, alt-forward and alt-reverse bases 計算每個變異位點的 AF
grep -v "^#" tmp.vcf | cut -f 8 | egrep -o 'DP4=[0-9]+,[0-9]+,[0-9]+,[0-9]+' | egrep -o '[0-9]+,[0-9]+,[0-9]+,[0-9]+' | awk -v FS="," '{print($3+$4)/($1+$2+$3+$4)}'
- 在前面步驟的bam文件里面找到這個vcf文件的某一個突變位點的測序深度表明的那些reads啼辣，并且在IGV里面可視化bam和vcf定位到該變異位點啊研。
less -SN tmp.vcf
截圖前幾行：
以34行和35行為例：
在BAM文件中定位reads:
samtools view tmp.bam | awk '{if ($6 != "*" && $4 <= 1104 && substr($10, 1104-$4+1, 1) == "A") print}' | less -SN
samtools view tmp.bam | awk '{if ($6 != "*" && $4 <= 1344 && substr($10, 1344-$4+1, 1) == "T") print}' | less -SN
IGV可視化：
其它思考題
- vcf的全稱是什么？是用來記錄什么信息的標(biāo)準(zhǔn)格式的文本？
Variant Call Format党远。記錄基因序列變異信息的文本文件削解。 - 一般選用哪個指令查看vcf文件，為什么不用vim?
less -SN沟娱，不用vim可能是處理大文件時效率低氛驮。 - vcf文件以’##’開頭的是什么信息？請認真查看這些信息济似〗梅希’#’開頭的是什么信息？
以 “##” 開頭的是注釋信息（meta-information）碱屁。以 “#” 開頭的是標(biāo)題行磷脯。 - Vcf文件除頭信息，每行有多少列娩脾，請詳細敘述每行的含義赵誓！請準(zhǔn)確記憶。
每行有8個固定的柿赊、必需的列俩功。第9列開始是基因型信息。
- 理解format列和樣本列的對應(yīng)關(guān)系以及GT AD DP的含義碰声。
format列和樣本列的數(shù)據(jù)是對應(yīng)的诡蜓，前者為格式，后者為格式對應(yīng)的數(shù)據(jù)胰挑。冒號分隔蔓罚。
GT: 樣本基因型。兩個數(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: Allele Depth，為sample中每一種allele的reads覆蓋度湃望，在二倍體中是用逗號分割的兩個值换衬，前者對應(yīng)ref基因型局义，后者對應(yīng)variant基因型。
DP: sample中該位點的覆蓋度冗疮。 - Vcf文件第三列如果不是’.’，出現(xiàn)的rs號的id是什么檩帐？
dbSNP數(shù)據(jù)庫中該SNP的ID號术幔。 - Vcf文件的ref，alt列和樣本列的0/1 1/1 或者1/2的聯(lián)系湃密？
- Vcf文件一般用什么軟件生成诅挑？請至少說出兩個軟件。請注意不同軟件生成的vcf格式的稍有不同的地方泛源。
bcftools拔妥，GATK4。 - Vcf文件一般都比較大达箍，壓縮后的.gz文件用什么指令直接查看而不用解壓后查看没龙？
zcat命令
gzip tmp.vcf
zcat tmp.vcf.gz
gunzip tmp.vcf.gz
- 了解gvcf是什么格式，gvcf全稱是什么缎玫？他與vcf有什么前后聯(lián)系硬纤？
Genome Variant Call Format。GVCF文件包含所有位點的記錄赃磨，既包括存在變異的位點筝家，也包括不存在變異的位點。這樣做的目的是為了方便后續(xù)的群體分析邻辉。 - 把alt列出現(xiàn)’,’的行提取出來
- 請將chrid溪王、postion、ref值骇、alt莹菱、format、樣本列切割出來生成一個文本
grep -v "^##" tmp.vcf | cut -f 1,2,4,5,9- > tmp_cut.txt
- 將一個含snp,indel信息的vcf拆成一個只含snp,一個只含indel信息的2個vcf文件雷客∶⒅椋可借鑒軟件
cat tmp.vcf | sed '/INDEL/d' > tmp_SNP.vcf # SNP
cat tmp.vcf | awk '/#|INDEL/{print $0}' > tmp_INDEL.vcf # INDEL
bcftools view tmp.vcf -O u -o tmp.bcf
bcftools view tmp.bcf
bcftools view -v snps tmp.bcf | less -SN # SNP
bcftools view -v indels tmp.bcf | less -SN # INDEL
vcftools --remove-indels --recode --recode-INFO-all --vcf tmp.vcf --stdout > tmp.snp.vcf # SNP
vcftools --keep-only-indels --recode --recode-INFO-all --vcf tmp.vcf --stdout > tmp.indel.vcf # INDEL
- 用指令操作indel的vcf文件，提取indel長度>4的變異行數(shù)搅裙，存成一個文本皱卓。
grep -v "^#" tmp.vcf | awk '{if (length($5) > 4) print}' > 4.txt
- 用Vcftools過濾vcf文件，如maf 設(shè)置成0.05部逮， depth設(shè)置成5-20娜汁，統(tǒng)計過濾前后的變異位點的總個數(shù)
vcftools --recode \
--recode-INFO-all \
--vcf tmp.vcf \
--maf 0.05 \
--minDP 5 --maxDP 20 \
--stdout > filter.vcf
# 會報錯，需要再回顧學(xué)習(xí)
- 利用vcftools提取每個樣本每一個位點的變異信息和深度信息兄朋，生成一個矩陣的文件掐禁，至少含義以下信息
Eg:
| CHRID | POSTION | SAMPLE_DP | SAMPLE_GT |
|---|---|---|---|
| chr1 | 1010 | 5 | AA |
- 提取出變異位點上樣本有純和突變的行數(shù)
- 統(tǒng)計一下1號染色體上的變異總個數(shù)。
- 提取一下BRCA1基因上發(fā)生變異的行數(shù)，如果是人類wes變異結(jié)果文件
- 統(tǒng)計vcf文件各樣本的缺失率傅事，如果是多個樣本的群體call結(jié)果缕允。
進階思考題：可視化vcf中變異位點質(zhì)量值，橫坐標(biāo)是質(zhì)量值蹭越，縱坐標(biāo)是百分比或者該質(zhì)量值時的變異位點的總個數(shù)障本，可以化成線圖或者點圖（可能需要學(xué)一些代碼，還有R畫圖）
參考：
http://www.bio-info-trainee.com/3577.html
https://samtools.github.io/hts-specs/VCFv4.2.pdf
http://www.reibang.com/p/eadfcd7e23ee
http://seqanswers.com/forums/showthread.php?t=10181
http://samtools.sourceforge.net/mpileup.shtml
https://cloud.tencent.com/developer/article/1054971
