linux test 比较两个数大小,linux_test

#1 在任意文件夹下面创建形如 1/2/3/4/5/6/7/8/9 格式的文件夹系列

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

#2 在创建好的文件夹下面，比如我的是 /Users/jimmy/tmp/1/2/3/4/5/6/7/8/9 ，里面创建文本文件 me.txt

#3 在文本文件 me.txt 里面输入内容:

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

I love bioinfomatics.

And you ?

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

#4 删除上面创建的文件夹 1/2/3/4/5/6/7/8/9 及文本文件 me.txt

rm -r 1/

#5 在任意文件夹下面创建 folder1~5这5个文件夹，然后每个文件夹下面继续创建 folder1~5这5个文件夹

mkdir -p folder{1..5}/folder{1..5}

#6 在第五题创建的每一个文件夹下面都 创建第二题文本文件 me.txt ，内容也要一样。

vim me.txt

echo folder{1..5}/folder{1..5} | xargs -n 1 cp me.txt

echo folder{1..5}/folder{1..5}/me.txt | xargs -n 1

图片.png

#7 再次删除掉前面几个步骤建立的文件夹及文件

rm -r folder{1..5}

#8 下载 http://www.biotrainee.com/jmzeng/igv/test.bed 文件，后在里面选择含有 H3K4me3 的那一行是第几行，该文件总共有几行。

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

grep -n H3K4me3 ~/test.bed

wc -l ~/test.bed

#9 下载 http://www.biotrainee.com/jmzeng/rmDuplicate.zip 文件，并且解压，查看里面的文件夹结构

#10 打开第九题解压的文件，进入 rmDuplicate/samtools/single 文件夹里面，查看后缀为 .sam 的文件，搞清楚 生物信息学里面的SAM/BAM 定义是什么。

less samtools/single/*.sam

#11 安装samtools

#12 打开 后缀为BAM 的文件，找到产生该文件的命令

samtools view -h tmp.sorted.bam | grep -n CL

#13 根据上面的命令，找到我使用的参考基因组 /home/jianmingzeng/reference/index/bowtie/hg38 具体有多少条染色体。

samtools view -h tmp.sorted.bam | grep SQ | awk '{print $2}' | cut -d"_" -f 1 | sort | uniq | wc -l

#14 上面的后缀为BAM 的文件的第二列，只有 0 和 16 两个数字，用 cut/sort/uniq等命令统计它们的个数。

samtools view -h tmp.sorted.bam | awk '{print $2}' | grep -v SN |sort -n | uniq -c

#15 重新打开 rmDuplicate/samtools/paired 文件夹下面的后缀为BAM 的文件，再次查看第二列，并且统计

samtools view -h tmp.sorted.bam | awk '{print $2}' | grep -v SN |sort -n | uniq -c

#16 下载 http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip 文件，并且解压，查看里面的文件夹结构， 这个文件有2.3M，注意留心下载时间及下载速度。

#17 解压 sickle-results/single_tmp_fastqc.zip 文件，并且进入解压后的文件夹，找到 fastqc_data.txt 文件，并且搜索该文本文件以 >>开头的有多少行？

cat fastqc_data.txt | grep ">>" | wc -l

#18 下载 http://www.biotrainee.com/jmzeng/tmp/hg38.tss 文件，去NCBI找到TP53/BRCA1等自己感兴趣的基因对应的 refseq数据库 ID，然后找到它们的hg38.tss 文件的哪一行。

cat hg38.tss | grep -n "NM_000546"

cat hg38.tss | grep -n NM_007294

#19 解析hg38.tss 文件，统计每条染色体的基因个数。

hg38.tss | awk '{print $2}' | cut -d"_" -f 1 | sort | uniq -c

#20 解析hg38.tss 文件，统计NM和NR开头的熟练，了解NM和NR开头的含义。

​cat hg38.tss | grep "NM" | wc -l

cat hg38.tss | grep "NR" | wc -l

图片.png

fasta&fastq

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

wc -l reads_1.fq

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

cat reads_1.fq | paste - - - - | awk '{print $1}' >1

or cat reads_1.fq | grep @r

or awk '{if(NR%4==1)print}' reads_1.fq >1 #(NR%4==1)每四行的第一行

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

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

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

awk '{if(NR%4==2)print}' reads_1.fq >2

awk '{if(NR%4==3)print}' reads_1.fq >3

awk '{if(NR%4==0)print}' reads_1.fq >4 #第四行是0！！!

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

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

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

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

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

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

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

https://www.jianshu.com/p/248308513e2e

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

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

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

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

awk '{if(NR%4==2)print}' reads_1.fq | cut -c1 | sort | uniq -c

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

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

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

wc -l reads1_1.fa

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

grep -o [GC] reads_1.fa | wc -l

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

cat reads_1.fa | tr -d 'N'

#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'

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

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

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

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

?

sam&bam

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

sed '1,3d' tmp.sam | wc -l

#2 比对类型看第二列

sed '1,3d' tmp.sam | awk '{print $2}' | sort -n | uniq -c

#3 比对成败看第三列

sed '1,3d' tmp.sam | awk '{if($3 == "*") print}'

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

single:

sed '1,3d' tmp.sam | awk '{if($3 == "*") print$1}' | sort | uniq -c |awk '{if($1 == "1") print$2}'| wc -l

sed '1,3d' tmp.sam | awk '{if($3 == "*") print$1}' | sort | uniq -c | grep -w 1| wc -l

double:

sed '1,3d' tmp.sam | awk '{if($3 == "*") print$1}' | sort | uniq -c |awk '{if($1 == "2") print$2}'| wc -l

sed '1,3d' tmp.sam | awk '{if($3 == "*") print$1}' | sort | uniq -c | grep -w 2| wc -l

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

sed '1,3d' tmp.sam | awk '$5 > 30 {print}'

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

sed '1,3d' tmp.sam | awk '{if($3 != "*") print}' | awk '{if($6 ~ /[IDNCSHP]/)print}' | wc -l

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

sed '1,3d' tmp.sam | awk '$9 > 1250 {print}'

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

sed '1,3d' tmp.sam | cut -f 3 | sort -u

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

sed '1,3d' tmp.sam | awk '{if($10 ~ /N/)print}'

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

sed '1,3d' tmp.sam | awk '{if($10 ~ /N/)print}' | awk '{if($3 != "*") print}' |

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

cat tmp.sam | grep '^@' | wc -l

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

cat tmp.sam | grep -v '^@' | cut -f 12- | awk -F" " '{print NF}' | sort | uniq -c

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

cat tmp.sam | grep -v '^@' | cut -f 12- | awk -F" " '{print NF}'

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

$grep "LN:" tmp.sam

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

sed '1,3d' tmp.sam | awk '{if($3 != "*") print}' | awk '{if($6 ~ /I/)print}' | less -S

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

sed '1,3d' tmp.sam | awk '{if($3 != "*") print}' | awk '{if($6 ~ /D/)print}' | less -S

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

cat tmp.sam |grep -v '^@'|awk '{if($4 > 5013 && $4 < 50130)print}'|less -S

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

cat tmp.sam | grep -v '^@'|sort -k4|less -S

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

grep 102M3D11M tmp.sam |cut -f 10|wc -m

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

概念题

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

SAM:The Sequence Alignment Map

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

可以直接查看

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

头部区：以’@'开始，体现了比对的一些总体信息。比如比对的SAM格式版本，比对的参考序列，比对使用的软件等。主体区：比对结果，每一个比对结果是一行，有11个主列和一个可选列。

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

以tab分割；有11个主列和一个可选列； awk -F ' ' '{print NF}'

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

序列ID及长度； grep @SQ tmp.sam | awk '{print $2 " " $3}'

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

从1开始

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.比对质量的数字是哪一列？越大比对质量越好还是越小越好？

第五列；越大越好

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

第二列；

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

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

samtools view -b ; bam文件为二进制，节省空间

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

samtools view; samtools view -h

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

BAM is compressed. Sorting helps to give a better compression ratio because similar sequences are grouped together.另外，samtools 对 BAM 文件进行排序之后那些没有比对上的 reads 会被放在文件的末尾。排序前头文件 SO:unsort, 排序后SO:coordinate； 排序后chr position是从1开始。

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

samtools index -b

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

IGV; 需要建立索引

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

samtools flagstat

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

单端测序得到的sam文件七八九列无意义(* 0 0)；

双端测序sam文件：第七列是这条reads第二次比对的位置，在利用bowtie2产生sam文件时，如果该列是“=”而第3列RNAME不是“*”则表示该reads两次比对均比对到第3列显示序列名的序列上；如果第3列RNAME和第7列MRNM都为“*”，则说明这条reads没有匹配上的序列；如果这条reads匹配到两个序列，则第一个序列的名称出现在第3列，而第二个序列的名称出现在第7列。

第8列表示与该reads对应的mate pair reads的比对位置，如果这对pair-end reads比对到同一条reference序列上，在sam文件中reads的id出现2次，Read1比对的第4列等于Read2比对的第8列。同样Read1比对的第8列等于Read2比对的第4列。

第9列可以理解为文库插入片段长度，如果R1端的read和R2端的read能够mapping到同一条Reference序列上(即第三列RNAME相同)，则该列的值表示第8列减去第4列加上第6列的值，R1端和R2端相同id的reads其第九列绝对值相同。

Bam能用gzip再压缩吗？

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

bwa; bowtie; boetie2

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

可以用sam2fasta或bam2fasta

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

bam文件直接打开后是二进制码

VCF

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

cat tmp.vcf | grep -v "^#" | grep INDEL | less -S

cat tmp.vcf | grep -v "^#" | grep -v INDEL | less -S

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

cat tmp.vcf | grep -v "^#" | grep INDEL | cut -f 8 | cut -d ";" -f 4 | cut -d"=" -f 2 |awk '{sum += $1};END {print sum/NR}'

cat tmp.vcf | grep -v "^#" | grep -v INDEL | cut -f 8 | cut -d ";" -f 4 | cut -d"=" -f 2 |awk '{sum += $1};END {print sum/NR}'

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

grep -v '^#' tmp.vcf | grep 'INDEL' | awk '{if(length($4) > length($5)) print }' | less -SN

grep -v '^#' tmp.vcf | grep 'INDEL' | awk '{if(length($4) < length($5)) print }' | less -SN

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

grep -v '^#' tmp.vcf | grep -v 'INDEL' | awk '{print $4$5}'| sort | uniq -c

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

qq | awk '{if ($10 !~ "1/1") print}' | wc -l

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

for i in `seq 21 100`; do grep -v '^#' tmp.vcf | cut -f8 | grep 'DP='${i} ; done | wc -l

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

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

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

for i in `seq 21 100`; do grep -v '^#' tmp.vcf | grep 'DP='${i} ; done | awk '{if ($6 > 30) print}' | wc -l

#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定位到该变异位点。

其他思考题

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

VCF是Variant Call Format的简称，是一种定义的专门用于存储基因序列突变信息的文本格式。例如基因组中的单碱基突变,SNP， 插入/缺失INDEL, 拷贝数变异CNV，和结构变异SV等，都是利用VCF格式来存储的。将其存储为二进制格式就是BCF。

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

less -S

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

##开头的是整体注释信息，#开头是各列的title

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

各列之间用tab空白隔开,前面9列为固定列，第10列开始为样品信息列，可以无限多个。

1.CHROM 记录染色体编号

2.POS 记录染色体位置信息

3.ID variant的ID。同时对应着dbSNP数据库中的ID，若没有，则默认使用‘.’; SNP/INDEL的dbSNP编号通常以rs开头，一般只有人类基因组才有dbSNP编号

4.REF 参考基因组碱基类型，必须是A,C,G,T,N且都大写。

5. ALT 变异碱基类型，必须是A,C,G,T,N且都大写，多个用逗号分割。"."表示这个地方没有reads覆盖为缺失。

6. QUAL 变异信息的检测质量值，越高越可靠。

7.FILTER 标记过滤结果的列，通常我们把VCF文件中的变异信息进行质控，过滤掉低质量的变异位点，如果该位点通过过滤标准那么我们可以在该列标记为"PASS",说明该列质量值高。标记完之后我们就可以用其他工具，把标记为"PASS"的列给筛选出来，这样方便后续分析。如果没有应用缺失值"."代替。

8.INFO 为附加信息列，一般以=；形式添加额外的注释信息列，常见的如DP=18 表示该位点测序深度为18X；AF=0.1表示等位基因频率为0.1

9.FORMAT 为后面10列信息的说明列，通常以":"隔开各个缩写词。

10.样品基因型列各信息以":"分隔与FORMAT列一一对应

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

GT 表示genotype，通常用”/” or “|”分隔两个数字；0代表参考基因组的碱基类型；1代表ALT碱基类型的第一个碱基(多个碱基用","分隔)，2代表ALT第二个碱基，以此类推；比如 REF列为：A， ALT列为G,T；那么0/1基因型为AG 杂合，1/1基因型为GG纯合SNP；1/2代表GT基因型；./.表示缺失；

AD: 两种碱基各自支持的碱基数量，用","分开两个数据，分别代表两个等位基因的深度；

DP:该样品该变异位点的测序深度总和，也就是AD两个数字的和；

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

SNP/INDEL的dbSNP编号

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

0代表参考基因组的碱基类型；1代表ALT碱基类型的第一个碱基(多个碱基用","分隔)，2代表ALT第二个碱基，以此类推；比如 REF列为：A， ALT列为G,T；那么0/1基因型为AG 杂合，1/1基因型为GG纯合SNP；1/2代表GT基因型；

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

bcftools， gatk

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

zcat

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

Genomic-Variant-Call-Format,gvcf文件与vcf文件都是vcf文件，不同之处在于gvcf文件会记录更多的信息，这里更多的信息指的是未突变的位点的覆盖情况

图片.png

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

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

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

grep -v "^##" tmp.vcf | awk '{print$2" "$3" "$4" "$5" "$9" "$10}' > test.txt

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

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

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

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

?

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

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