20181204

作者: 0508aa33893c | 来源:发表于2018-12-04 16:21 被阅读108次

########################## 生信人的linux考试 ##########################

#### 1 创建层级文件夹

mkdir -p ~/1/2/3/4/5/6/7/8/9 && cd ~/1/2/3/4/5/6/7/8/9

#### 2 指定地方创建文本文件

vim me.txt

#### 3 写入文本文件

按I键，输入如下内容：

Go to: http://www.biotrainee.com/

I love bioinfomatics.

And you ?

输入完毕后，按ESC ：wq

cat me.txt

cp ~/1/2/3/4/5/6/7/8/9/me.txt ~/

#### 4

rm -r ~/1

#### 5

mkdir -pv folder_{1..5}/folder_{1..5}

tree

cd ~

#### 6

vim cp.sh

#!/bin/bash

for I in {1..5};do

cp ~/me.txt ~/folder_$I/

for i in {1..5};do

cp ~/me.txt ~/folder_$I/folder_$i

done

chmod u+x cp.sh

~/cp.sh

#### 7

rm -r folder_*

#### 8 找出

wget -c http://www.biotrainee.com/jmzeng/igv/test.bed

grep -n -o --color H3K4me3 test.bed ==> 第8行

wc -l test.bed ==> 共10行

#### 9

wget -c http://www.biotrainee.com/jmzeng/rmDuplicate.zip

unzip rmDuplicate.zip && ls -l

tree rmDuplicate

#### 10

cd ~/rmDuplicate/samtools/single

less -SN tmp.sam

cat tmp.sam | head

#sam 是序列比对的标准格式,头部加比对，比对12列以tab分隔，12列包括：1序列名称 2标记 3比对上的参考序列的染色体号 4染色体起点位置 5质量分数 6CRGAI的字符串换 789：单端测序789行没有 7mate比对到的染色体号 8mate比对到参考序列上的第一个碱基位置 9两条read之间的间隔 10 read的序列信息 11 ASCII序列质量 12 可选字段

#bam是二进制，占内存要小，节省存储空间

#### 11

cd ~/biosoft

wget -c https://sourceforge.net/projects/samtools/files/samtools/1.9/samtools-1.9.tar.bz2

tar -jxvf samtools-1.9.tar.bz2

cd samtools-1.9

make

#### 12

cd ~/rmDuplicate/samtools/single

vip48 23:49:26 ~/linux_test_20/rmDuplicate/samtools/single

vim ~/.bashrc

export PATH="$PATH:~/biosoft/samtools/samtools-1.9"

source ~/.bashrc

samtools view tmp.rmdup.bam

cat tmp.header

找到命令CL:"/home/jianmingzeng/biosoft/bowtie/bowtie2-2.2.9/bowtie2-align-s --wrapper basic-0 -p 20 -x /home/jianmingzeng/reference/index/bowtie/hg38 -S /home/jianmingzeng/data/public/allMouse/alignment/WT_rep2_Input.sam -U /tmp/41440.unp"

#### 13

less tmp.header | head

#### 14

samtools view -H sorted.bam|awk '{print $2}'|cut -c4-9|grep -v '_' |sort|uniq -cd

samtools view tmp.rmdup.bam | cut -f2 | sort -n |uniq -dc ==> 0有16个，16有12个

#### 15

cd ~/rmDuplicate/samtools/paired

samtools view tmp.rmdup.bam | cut -f2 | sort -n |uniq -dc ==> 83有2个，97有2个，99有8个，147有7个，163有2个

#### 16

wget -c http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip ==>速度80 kb/s,

unzip sickle-results.zip

tree

#### 17

cd sickle-results/

unzip single_tmp_fastqc.zip

cd single_tmp_fastqc/

grep ^\>\> fastqc_data.txt | wc -l ==> 24行

#### 18

wget -c http://www.biotrainee.com/jmzeng/tmp/hg38.tss

grep NM_001126113 hg38.tss ==> 29346行

#### 19

cat hg38.tss | cut -f2 | cut -c1-5 | sort | uniq -dc

#### 20

grep -o -E "NM|NR" hg38.tss |sort |uniq -dc

NM:转录组产物的序列 ==> 51064

NR:非编码的转录组序列==> 15954

################## fasta和fastq格式文件的shell小练习 ##################

mkdir -p ~/biosoft

cd ~/biosoft

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

1) 统计reads_1.fq 文件中共有多少条序列信息

cat reads_1.fq | paste - - - - | wc -l ==> 10000条序列

2) 输出所有的reads_1.fq文件中的标识符(即以@开头的那一行)

cat reads_1.fq | paste - - - - | cut -f1 | tee reads_1.f1.fq

3) 输出reads_1.fq文件中的所有序列信息(即每个序列的第二行)

cat reads_1.fq | paste - - - - | cut -f2 | tee reads_1.f2.fq

4) 输出以‘+’及其后面的描述信息(即每个序列的第三行)

cat reads_1.fq | paste - - - - | cut -f3 | tee reads_1.f3.fq

5) 输出质量值信息(即每个序列的第四行)

cat reads_1.fq | paste - - - - | cut -f4 | tee reads_1.f4.fq

6) 计算reads_1.fq 文件含有N碱基的reads个数

cat reads_1.fq | paste - - - - | cut -f2 | grep -c "N" ==> 6429个reads # 行数目

7) 统计文件中reads_1.fq文件里面的序列的碱基总数

cat reads_1.fq | paste - - - - | cut -f2 | grep -E -o "A|T|C|G|N" | wc -l ==> 1088399个碱基

awk '{if(NR%4==2)print}' reads_1.fq |grep -o [ATCGN] | wc -l

awk '{if(NR%4==2)print length}' reads_1.fq | paste -s -d + | bc

8) 计算reads_1.fq 所有的reads中N碱基的总数

cat reads_1.fq | paste - - - - | cut -f2 | grep -o "N" | wc -l ==> 26001 个

9) 统计reads_1.fq 中测序碱基质量值恰好为Q20的个数

cat reads_1.fq | paste - - - - | cut -f4 | grep -o "5" | wc -l ==> 21369

10) 统计reads_1.fq 中测序碱基质量值恰好为Q30的个数

cat reads_1.fq | paste - - - - | cut -f4 | grep -o "?" | wc -l ==> 21574

11) 统计reads_1.fq 中所有序列的第一位碱基的ATCGNatcg分布情况

==> 2184 A

==> 2203 C

==> 2219 G

==> 1141 N

==> 2253 T

12) 将reads_1.fq 转为reads_1.fa文件(即将fastq转化为fasta)

$cat reads_1.fq | paste - - - - | cut -f1-2 |tr '@' '>' | tr '\t' '\n' >> reads_1.fa

13) 统计上述reads_1.fa文件中共有多少条序列

cat reads_1.fa | paste - - | wc -l ==> 10000

14) 计算reads_1.fa文件中总的碱基序列的GC数量

cat reads_1.fa | paste - - | cut -f2 | grep -o [GC] | wc -l ==> 529983

15) 删除 reads_1.fa文件中的每条序列的N碱基

cat reads_1.fa | paste - - | cut -f2 | tr -d 'N' >> reads_1.del.fa

16) 删除 reads_1.fa文件中的含有N碱基的序列

cat reads_1.fa | paste - - | grep -v "N" | tr '\t' '\n'

17) 删除 reads_1.fa文件中的短于65bp的序列

cat reads_1.fa | paste - - | awk '{ if(length($2)>=65)print }' | tr '\t' '\n' >> reads_1.65.fa

18) 删除 reads_1.fa文件每条序列的前后五个碱基

cat reads_1.fa | paste - - |awk '{$a=substr($2,6,length($2)-10);print $1,$a;}'OFS='\t' >> reads_1.5.fa

19) 删除 reads_1.fa文件中的长于125bp的序列

cat reads_1.fa | paste - - | awk '{ if(length($2)<=125)print }' | tr '\t' '\n' >> reads_1.125.fa

20) 查看reads_1.fq 中每条序列的第一位碱基的质量值的平均值

#cat reads_1.fq | paste - - - - | cut -f4 |cut -c1|tr "[!"#$%&,()*+,-./0123456789:;<=>?@ABCDEFGHI]" "[33-73]"|awk '{sum+=$1} END {print sum/NR}'|tr "[33-73]" "[!"#$%&,()*+,-./0123456789:;<=>?@ABCDEFGHI]"

？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？

？？？？？？？？

？？？？？？？？？？？？？？？？？？

################## sam和bam格式文件的shell小练习 ##################

cd ~/biosoft

wget -c 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 > tmp.sam

samtools view -bS tmp.sam >tmp.bam

1) 统计共多少条reads(pair-end reads这里算一条)参与了比对参考基因组

samtools view tmp.sam | cut -f1 |sort|uniq|wc -l ==> 10000条reads

2) 统计共有多少种比对的类型(即第二列数值有多少种)及其分布

samtools view tmp.sam | cut -f2 |sort|uniq -dc|sort -k1 -nr

==> 结果如下：

4650 83

4650 163

4516 99

4516 147

213 77

213 141

165 69

165 137

153 73

153 133

136 89

136 165

125 153

125 101

24 65

24 129

16 177

16 113

2 81

2 161

3）筛选出比对失败的reads，看看序列特征

samtools view tmp.sam | awk '{if($6=="*")print $10}' | less -SN ==> 1005个好多N

4) 比对失败的reads区分成单端失败和双端失败情况，并且拿到序列ID

单端：

samtools view tmp.sam | awk '{if($6=="*")print $1}'|sort|uniq -c|grep -w 1

双端：

samtools view tmp.sam | awk '{if($6=="*")print $1}'|sort|uniq -c|grep -w 2

5) 筛选出比对质量值大于30的情况（看第5列）

samtools view tmp.sam | awk '{if($5>30)print}' > tmp.50.sam

6) 筛选出比对成功，但是并不是完全匹配的序列

samtools view tmp.sam | awk '{if($6!="*" && $6 ~"[IDNXSHP]")print}'

#samtools view tmp.sam | awk '{if($6!="*")print $6}'|grep "[IDNXSHP]"

7) 筛选出inset size长度大于1250bp的 pair-end reads

samtools view tmp.sam | awk '{if($9>1250 || $9<-1250)print $0}'|less -S

8) 统计参考基因组上面各条染色体的成功比对reads数量

samtools view tmp.sam | cut -f3 | sort|uniq -c ==> 19574 gi|9626243|ref|NC_001416.1|

9) 筛选出原始fq序列里面有N的比对情况

samtools view tmp.sam |awk '{if($10 ~"N")print}'

10) 筛选出原始fq序列里面有N，但是比对的时候却是完全匹配的情况

samtools view tmp.sam |awk '{if($10 ~"N" && $6 !~"[IDNXSHP*]")print}'

11) sam文件里面的头文件行数

cat tmp.sam | grep "^@" | wc -l ==> 3行

12) sam文件里每一行的tags个数一样吗

# samtools view tmp.sam|head -1|tr "\t" "\n"|wc -l

samtools view tmp.sam| awk '{print NF-11}'|sort|uniq -c ==> 不一样,如下：

457 1

548 2

579 8

18416 9

13) sam文件里每一行的tags个数分别是多少个

samtools view tmp.sam| awk '{print $1"\t"NF-11}'

14) sam文件里记录的参考基因组染色体长度分别是？

samtools view -H tmp.sam |grep "LN:" ==> 48502

15) 找到比对情况有insertion情况的

samtools view tmp.sam |awk '{if($6 ~"I")print}'

16) 找到比对情况有deletion情况的

samtools view tmp.sam |awk '{if($6 ~"D")print}'

17）取出位于参考基因组某区域的比对记录，比如 5013到50130 区域

samtools view tmp.sam |awk '{if($4<50130 && $4>5013)print}'

18) 把sam文件按照染色体以及起始坐标排序

samtools view tmp.sam |sort -k3,4

19) 找到 102M3D11M 的比对情况，计算其reads片段长度

samtools view tmp.sam |grep "102M3D11M" | awk '{print length($10)}' ==> 113

20) 安装samtools软件后使用samtools软件的各个功能尝试把上述题目重新做一遍

？？？？？？？？？？？？还没学软件安装，接下来要看

其它概念题

1）高通量测序数据分析中，序列与参考序列的比对后得到的标准文件为什么文件？

==> sam文件

2）sam文件是一种比对后的文件，能直接查看吗？

==> 可以

3）Sam/Bam文件分为两部分，一部分以@开头的是什么序列，另一部分是什么序列？

==> @开头的是注释信息头文件另一部分是比对结果部分

4）Sam文件除头文件，以什么符号分割文本的，比对部分信息一行是多少列？你能用awk计算列数吗？

==> 以tab制表符分割文件，比对部分信息一行除了可选字段有11列，awk '{print NF}'

5）Sam/Bam文件的@SQ开头的行是什么？你能生成一个文本，两列，一列是参考基因组染色体/sca id，一列是长度吗？

==> 是reads比对上的染色体及染色体的长度，cat tmp.sam|grep -w ^@SQ|cut -f2-3 > chrom.txt

6）Sam文件的比对位置是从1还是0开始的？

==> 从1开始，0表示没比对上

7）常见CIGAR 字符串各字母代表的意义

==>

M：alignment match (can be a sequence match or mismatch)

表示read可mapping到第三列的序列上，则read的碱基序列与第三列的序列碱基相同，表示正常的mapping结果，M表示完全匹配，但是无论reads与序列的正确匹配或是错误匹配该位置都显示为M

I：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”表示错误匹配到序列上

8）比对质量的数字是哪一列？越大比对质量越好还是越小越好？

==> 第5列，数字越大特异性越高越好

9）Sam文件的flag是第几列，数字代表什么意义？是怎么计算来的？

==> 第2列，

数值结果如下：

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是次优的比对结果

512（1000000000）该read没有通过质量控制

1024（10000000000）由于PCR或测序错误产生的重复reads

2048（100000000000）补充匹配的read

具体的flag值的解释，可以参考samtools软件提供的结果

符合特质的数值的加和得来

10）Sam文件怎么转bam文件？用什么指令？为什么要转换？

==> samtools view -bS tmp.sam >tmp.bam bam占用的磁盘空间比sam文本文件小；利用bam二进制文件的运算速度快

11）Bam文件查看用什么指令？如果需要查看头文件需要增加参数？

==> samtools view tmp.bam samtools -H view tmp.bam

12）Bam文件为什么要排序？排序后的bam和未排序的bam头文件和chr position列有什么区别？

==> 排序后会按照染色体位置排序排好

13）Bam文件建索引的指令是什么？

==> samtools index <in.bam> [out.index]

14）Bam文件可视化用什么工具？查看时需要建立索引吗？

==> IGV 需要

15）Bam文件统计reads比对情况用samtools的哪一个子命令？

==> samtools idxstats <aln.bam>

16）SE测序和PE测序的所比对后得到的sam文件的区别在哪里？

==> Qname第一列，双端有重复现象,单端没有重复

17）Bam能用gzip再压缩吗？

==> 不能，tmp.bam.gz already has .gz suffix -- unchanged

18）Sam文件通常由哪些比对软件得到的

==> Bowtie2、BWA、

19）Sam/Bam文件能转成fastq文件吗？

==> 能

samtools view -bS aln.sam > aln.bam

samtools view -h -o aln.sam aln.bam

bam2fastq --aligned input.bam -o output.fq

20）有时候不能通过文件名的后缀来区别是sam还是bam，你是怎么区别的？

==> cat后，bam是乱码

################## VCF格式文件的shell小练习 ##################

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 -c 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练习题

1) 把突变记录的vcf文件区分成 INDEL和SNP条目#

cat tmp.vcf|grep -v "^#"|grep INDEL ==> INDEL

cat tmp.vcf|grep -v "^#"|grep -v INDEL ==> SNP

2) 统计INDEL和SNP条目的各自的平均测序深度

cat tmp.vcf|grep -v "^#"|grep INDEL|cut -f8|cut -d";" -f4|cut -d"=" -f2|awk '{ sum += $0; } END { print sum/NR }' ==> 40.4429 <== INDEL

cat tmp.vcf|grep -v "^#"|grep -v INDEL|cut -f8|cut -d";" -f1|cut -d"=" -f2|awk '{ sum += $0; } END { print sum/NR }' ==> 37 <== SNP

3) 把INDEL条目再区分成insertion和deletion情况

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

4) 统计SNP条目的突变组合分布频率

cat tmp.vcf|grep -v "^#"|grep -v INDEL|awk '{print $4"-->"$5}'|sort|uniq -c|sort -k1,1nr

5) 找到基因型不是 1/1 的条目，个数

less tmp.vcf|grep -v "^#"|grep 1/1

less tmp.vcf|grep -v "^#"|grep 1/1|wc -l ==> 81个

6) 筛选测序深度大于20的条目

cat tmp.vcf|grep -v "^#"|awk '{if($DP>20)print}'

7) 筛选变异位点质量值大于30的条目

cat tmp.vcf|grep -v "^#"|awk '{if($6>30)print}'

8) 组合筛选变异位点质量值大于30并且深度大于20的条目

cat tmp.vcf|grep -v "^#"|awk '{if($6>30 && $DP>20)print}'

9) 理解DP4=4,7,11,18 这样的字段，就是 Number of high-quality ref-forward , ref-reverse, alt-forward and alt-reverse bases 计算每个变异位点的 AF

???????????????????? ???????????

10）在前面步骤的bam文件里面找到这个vcf文件的某一个突变位点的测序深度表明的那些reads，并且在IGV里面可视化bam和vcf定位到该变异位点。

??????????????????????? IGV ?????????

其它思考题

1）vcf的全称是什么？是用来记录什么信息的标准格式的文本？

VCF文件全称为Variant Call Format，表示基因组的变异信息，通常为GATK和Samtools软件处理所得到，vcf是用于描述SNP,INDEL和SV的结果的文本文件，VCF全称为Variant Call Format, 是一种纯文本文件，用来存储变异位点信息，主要包括3个部分的内容

## mate-information line，# header line，data line

2）一般选用哪个指令查看vcf文件，为什么不用vim?

???????????????????????

3）vcf文件以’##’开头的是什么信息？请认真查看这些信息。’#’开头的是什么信息？

## mate-information line，以##开头，格式为key=value。 fileformat是必须的字段，表明VCF格式的版本，其他行主要用来描述INFO, FORMAT, FILTER等字段的具体含义。

# header line,以#开头，\t分隔，至少包含8个字段

4）Vcf文件除头信息，每行有多少列，请详细叙述每行的含义！请准确记忆。

至少9列

CHROM : 染色体名字

POS : 染色体的位置，起始位置为1

ID : 变异位点在数据库中的ID,如果是dbsnp数据库，推荐使用rs号，如果没有ID,用点号表示缺失值

REF : 参考基因组上的碱基

ALT : 变异之后的碱基

QUAL : 变异位点的质量，质量值越高，为真实的变异位点的概率越大

FILTER : 过滤信息，PASS 代表通过了过滤；对于过滤失败的位点，会给出对应的过滤失败的原因，具体的含义可以查看mate-information line 中对FILTER 字段的描述

INFO ：额外的信息，具体的含义可以查看mate-information line 中对INFO 字段的描述,看上去是一列，但其中的内容可以无限扩增

Format:表头的##FORMAT就是对第九列的解释，主要包括某一个特定位点基因型、测序深度的表述

5）理解format列和样本列的对应关系以及GT AD DP的含义。

GT：样本基因型（genotype），两个数字之间【这里是1/1】斜线分隔，表示二倍体样本的基因型。

AD和DP：第八列也有DP，但含义不同

AD是Allele Depth，样本中每一种allel的reads覆盖度，在二倍体中是用逗号分隔的两个数，前面对应ref，后面对应variant；

DP是Depth，是样本中该位点覆盖度

6）Vcf文件第三列如果不是’.’，出现的rs号的id是什么？

ID：变异位点名称,对应dbSNP数据库中的ID

7）Vcf文件的ref，alt列和样本列的0/1 1/1 或者1/2的联系？

0代表样本中ref的allel，

1代表样本variant的allel，

2表示有第二个variant的allel

0/0表示样本中该位点纯合，与ref一致；

0/1表示样本中该位点杂合，有ref和variant两个基因型；

1/1表示样本中位点纯合，与variant一致

8）Vcf文件一般用什么软件生成？请至少说出两个软件。请注意不同软件生成的vcf格式的稍有不同的地方。

通常为GATK和Samtools软件处理所得到

9）Vcf文件一般都比较大，压缩后的.gz文件用什么指令直接查看而不用解压后查看？

zless zcat

10）了解gvcf是什么格式，gvcf全称是什么？他与vcf有什么前后联系？

GVCF代表Genomic VCF,GVCF是一种VCF，因此基本格式规范与常规VCF相同，但基因组VCF包含额外信息。

术语GVCF有时仅用于描述包含基因组中每个位置（或感兴趣区间）的记录的VCF，

无论是否在该位点检测到变体（例如由UnifiedGenotyper产生的VCF --output_mode EMIT_ALL_SITES）。

HaplotypeCaller 3.x生成的GVCF包含以非常特定的方式格式化的其他信息。

?????????????????????

11）把alt列出现’,’的行提取出来

cat tmp.vcf |grep -v "^#"|awk '{if($5 ~",")print}'

12）请将chrid、postion、ref、alt、format、样本列切割出来生成一个文本

cat tmp.vcf |grep -v "^#"|awk '{print $1,$2,$4,$5,$9}'> tmp.new.txt

13）将一个含snp,indel信息的vcf拆成一个只含snp,一个只含indel信息的2个vcf文件。可借鉴软件

?????????????????? 还不会用软件 😭

14）用指令操作indel的vcf文件，提取indel长度>4的变异行数，存成一个文本。

15）用Vcftools过滤vcf文件，如maf 设置成0.05， depth设置成5-20，统计过滤前后的变异位点的总个数

16）利用vcftools提取每个样本每一个位点的变异信息和深度信息，生成一个矩阵的文件，至少含义以下信息

Eg:

CHRID POSTION SAMPLE_DP SAMPLE_GT

chr1 1010 5 AA

17）提取出变异位点上样本有纯和突变的行数

18）统计一下1号染色体上的变异总个数。

19）提取一下BRCA1基因上发生变异的行数，如果是人类wes变异结果文件

20）统计vcf文件各样本的缺失率，如果是多个样本的群体call结果。

21）进阶思考题：可视化vcf中变异位点质量值，横坐标是质量值，纵坐标是百分比或者该质量值时的变异位点的总个数，可以化成线图或者点图（可能需要学一些代码，还有R画图）

20181204

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

Linux基础