首先使用bowtie2软件自带的测试数据生成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练习题
- 把突变记录的vcf文件区分成 INDEL和SNP条目
grep 'INDEL' tmp.vcf
grep -v 'INDEL' tmp.vcf
- 统计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}'
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}'
- 把INDEL条目再区分成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
- 统计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 的条目，个数
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 | wc -l
- 筛选变异位点质量值大于30的条目
grep -v '^#' tmp.vcf | awk '{if ($6 > 30) print}' | less -SN
- 组合筛选变异位点质量值大于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 | awk '{print $10}' | awk -F ":" '{print$1}' | sort | uniq -c
grep -v "^#" tmp.vcf | awk '{if ($0 ~ "1/1") {$8=$8";AF=1"; print$0} else if ($0 ~ "0/1") {$8=$8";AF=0.5"; print$0}}' | less -SN
- 在前面步骤的bam文件里面找到这个vcf文件的某一个突变位点的测序深度表明的那些reads，并且在IGV里面可视化bam和vcf定位到该变异位点。
其它思考题
vcf的全称是什么？是用来记录什么信息的标准格式的文本？
一般选用哪个指令查看vcf文件，为什么不用vim?
vcf文件以’##’开头的是什么信息？请认真查看这些信息。’#’开头的是什么信息？
Vcf文件除头信息，每行有多少列，请详细叙述每行的含义！请准确记忆。
理解format列和样本列的对应关系以及GT AD DP的含义。
Vcf文件第三列如果不是’.’，出现的rs号的id是什么？
Vcf文件的ref，alt列和样本列的0/1 1/1 或者1/2的联系？
Vcf文件一般用什么软件生成？请至少说出两个软件。请注意不同软件生成的vcf格式的稍有不同的地方。
Vcf文件一般都比较大，压缩后的.gz文件用什么指令直接查看而不用解压后查看？
了解gvcf是什么格式，gvcf全称是什么？他与vcf有什么前后联系？
把alt列出现’,’的行提取出来
请将chrid、postion、ref、alt、format、样本列切割出来生成一个文本
将一个含snp,indel信息的vcf拆成一个只含snp,一个只含indel信息的2个vcf文件。可借鉴软件
用指令操作indel的vcf文件，提取indel长度>4的变异行数，存成一个文本。
用Vcftools过滤vcf文件，如maf 设置成0.05， depth设置成5-20，统计过滤前后的变异位点的总个数
利用vcftools提取每个样本每一个位点的变异信息和深度信息，生成一个矩阵的文件，至少含义以下信息
Eg:
| CHRID | POSTION | SAMPLE_DP | SAMPLE_GT |
|---|---|---|---|
| chr1 | 1010 | 5 | AA |
提取出变异位点上样本有纯和突变的行数
统计一下1号染色体上的变异总个数。
提取一下BRCA1基因上发生变异的行数，如果是人类wes变异结果文件
统计vcf文件各样本的缺失率，如果是多个样本的群体call结果。
进阶思考题：可视化vcf中变异位点质量值，横坐标是质量值，纵坐标是百分比或者该质量值时的变异位点的总个数，可以化成线图或者点图（可能需要学一些代码，还有R画图）
参考：
http://www.bio-info-trainee.com/3577.html
https://www.jianshu.com/p/eadfcd7e23ee
网友评论