首先使用bowtie2軟件自帶的測(cè)試數(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í)題

1. 把突變記錄的vcf文件區(qū)分成 INDEL和SNP條目
2. 統(tǒng)計(jì)INDEL和SNP條目的各自的平均測(cè)序深度
3. 把INDEL條目再區(qū)分成insertion和deletion情況
4. 統(tǒng)計(jì)SNP條目的突變組合分布頻率
5. 找到基因型不是1/1的條目，個(gè)數(shù)
6. 篩選測(cè)序深度大于20的條目
7. 篩選變異位點(diǎn)質(zhì)量值大于30的條目
8. 組合篩選變異位點(diǎn)質(zhì)量值大于30并且深度大于20的條目
9. 理解DP4=4,7,11,18這樣的字段，就是 Number of high-quality ref-forward , ref-reverse, alt-forward and alt-reverse bases 計(jì)算每個(gè)變異位點(diǎn)的 AF
10. 在前面步驟的bam文件里面找到這個(gè)vcf文件的某一個(gè)突變位點(diǎn)的測(cè)序深度表明的那些reads，并且在IGV里面可視化bam和vcf定位到該變異位點(diǎn)。

其它思考題

1. vcf的全稱是什么？是用來(lái)記錄什么信息的標(biāo)準(zhǔn)格式的文本？

2. 一般選用哪個(gè)指令查看vcf文件，為什么不用vim?

3. vcf文件以’##’開(kāi)頭的是什么信息？請(qǐng)認(rèn)真查看這些信息。’#’開(kāi)頭的是什么信息？

4. Vcf文件除頭信息，每行有多少列，請(qǐng)?jiān)敿?xì)敘述每行的含義！請(qǐng)準(zhǔn)確記憶。

5. 理解format列和樣本列的對(duì)應(yīng)關(guān)系以及GT AD DP的含義。

6. Vcf文件第三列如果不是’.’，出現(xiàn)的rs號(hào)的id是什么？

7. Vcf文件的ref，alt列和樣本列的0/1 1/1 或者1/2的聯(lián)系？

8. Vcf文件一般用什么軟件生成？請(qǐng)至少說(shuō)出兩個(gè)軟件。請(qǐng)注意不同軟件生成的vcf格式的稍有不同的地方。

9. Vcf文件一般都比較大，壓縮后的.gz文件用什么指令直接查看而不用解壓后查看？

10. 了解gvcf是什么格式，gvcf全稱是什么？他與vcf有什么前后聯(lián)系？

11. 把a(bǔ)lt列出現(xiàn)’,’的行提取出來(lái)

12. 請(qǐng)將chrid、postion、ref、alt、format、樣本列切割出來(lái)生成一個(gè)文本

13. 將一個(gè)含snp,indel信息的vcf拆成一個(gè)只含snp,一個(gè)只含indel信息的2個(gè)vcf文件?？山梃b軟件

14. 用指令操作indel的vcf文件，提取indel長(zhǎng)度>4的變異行數(shù)，存成一個(gè)文本。

15. 用Vcftools過(guò)濾vcf文件，如maf 設(shè)置成0.05， depth設(shè)置成5-20，統(tǒng)計(jì)過(guò)濾前后的變異位點(diǎn)的總個(gè)數(shù)

16. 利用vcftools提取每個(gè)樣本每一個(gè)位點(diǎn)的變異信息和深度信息，生成一個(gè)矩陣的文件，至少含義以下信息

Eg:

1. 提取出變異位點(diǎn)上樣本有純和突變的行數(shù)

2. 統(tǒng)計(jì)一下1號(hào)染色體上的變異總個(gè)數(shù)。

3. 提取一下BRCA1基因上發(fā)生變異的行數(shù)，如果是人類wes變異結(jié)果文件

4. 統(tǒng)計(jì)vcf文件各樣本的缺失率，如果是多個(gè)樣本的群體call結(jié)果。

進(jìn)階思考題：可視化vcf中變異位點(diǎn)質(zhì)量值，橫坐標(biāo)是質(zhì)量值，縱坐標(biāo)是百分比或者該質(zhì)量值時(shí)的變異位點(diǎn)的總個(gè)數(shù)，可以化成線圖或者點(diǎn)圖（可能需要學(xué)一些代碼，還有R畫(huà)圖）

本文作者：生信技能樹(shù)團(tuán)隊(duì)