LINUX练习题
1) 统计共多少条reads(pair-end reads这里算一条)参与了比对参考基因组
$grep -v '^@' tmp.sam |cut -f 1|sort -V | uniq -c | wc -l
10000
- 统计共有多少种比对的类型(即第二列数值有多少种)及其分布。
$grep -v '^@' tmp.sam | cut -f 2 | sort -V |uniq -c
24 65
165 69
153 73
213 77
2 81
4650 83
136 89
4516 99
125 101
16 113
24 129
153 133
165 137
213 141
4516 147
125 153
2 161
4650 163
136 165
16 177
$grep -v '^@' tmp.sam |cut -f 2|sort -V | uniq -c | wc -l
20
3）筛选出比对失败的reads，看看序列特征。
$grep -v '^@' tmp.sam | awk 'match($3,"*"){print;}'
$grep -v '^@' tmp.sam | awk 'match($3,"*"){print;}' | wc -l
426
$grep -v '^@' tmp.sam | awk 'match($3,"*"){print;}' | head -10 |less -SN
- 比对失败的reads区分成单端失败和双端失败情况，并且拿到序列ID
$grep -v '^@' tmp.sam | awk 'match($3,"*"){print($1);}' |\
sort -V | uniq -c | awk '$1==2{print($2);}' #双端失败
$grep -v '^@' tmp.sam | awk 'match($3,"*"){print($1);}' |\
sort -V | uniq -c | awk '$1==2{print($2);}' | wc -l
213
$grep -v '^@' tmp.sam | awk 'match($3,"*"){print($1);}' |\
sort -V | uniq -c | awk '$1==1{print($2);}' #单端失败
$grep -v '^@' tmp.sam | awk 'match($3,"*"){print($1);}' |\
sort -V | uniq -c | awk '$1==1{print($2);}' | wc -l
0
- 筛选出比对质量值大于30的情况（看第5列）
$grep -v '^@' tmp.sam | awk '$5>30{print;}'
$grep -v '^@' tmp.sam | awk '$5>30{print;}' | wc -l
18632
$grep -v '^@' tmp.sam | awk '$5>30{print;}'| head -10 |less -SN
- 筛选出比对成功，但是并不是完全匹配的序列
- 筛选出inset size长度大于1250bp的 pair-end reads
- 统计参考基因组上面各条染色体的成功比对reads数量
$grep -v '^@' tmp.sam | awk '!match($6,"*"){print($3);}' | sort | uniq -c
18995 gi|9626243|ref|NC_001416.1|
- 筛选出原始fq序列里面有N的比对情况
$grep -v '^@' tmp.sam | awk 'match($10,"N"){print}'
$grep -v '^@' tmp.sam | awk 'match($10,"N"){print}'| wc -l
12934
$grep -v '^@' tmp.sam | awk 'match($10,"N"){print}' | \
head -10 | cut -f 10 | grep "N"
- 筛选出原始fq序列里面有N，但是比对的时候却是完全匹配的情况
$grep -v '^@' tmp.sam | awk 'match($6,length($10)){print}'
$grep -v '^@' tmp.sam | awk 'match($6,length($10)){print}'| wc -l
17095
$grep -v '^@' tmp.sam | awk 'match($6,length($10)){print}'| head -10 | less -SN
- sam文件里面的头文件行数
$grep '^@' tmp.sam | wc -l
3
- sam文件里每一行的tags个数一样吗
$grep -v '^@' tmp.sam |awk '{print(NF-12+1)}'| sort | uniq -c
457 1
548 2
579 8
18416 9
- sam文件里每一行的tags个数分别是多少个
$grep -v '^@' tmp.sam |awk '{print(NF-12+1)}'
$grep -v '^@' tmp.sam |awk '{print(NF-12+1)}'|head -1
9
$grep -v '^@' tmp.sam |head -1 |\
awk '{for(i=12;i<NF;i++){printf("%s\n",$i)}}' | less -SN
013-2.png
- sam文件里记录的参考基因组染色体长度分别是？
$grep "LN:" tmp.sam
@SQ SN:gi|9626243|ref|NC_001416.1| LN:48502
- 找到比对情况有insertion情况的
$grep -v '^@' tmp.sam | awk '{if (match($6,"[Ii]")) print;}'
$grep -v '^@' tmp.sam | awk '{if (match($6,"[Ii]")) print;}'| wc -l
66
$grep -v '^@' tmp.sam | awk '{if (match($6,"[Ii]")) print;}'|\
head -10 | less -SN
015.png
- 找到比对情况有deletion情况的
$grep -v '^@' tmp.sam | awk '{if (match($6,"[Dd]")) print;}'
$grep -v '^@' tmp.sam | awk '{if (match($6,"[Dd]")) print;}'| wc -l
1843
$grep -v '^@' tmp.sam | awk '{if (match($6,"[Dd]")) print;}'|\
head -10 | less -SN
17）取出位于参考基因组某区域的比对记录，比如 5013到50130 区域
$grep -v '^@' tmp.sam | awk '{if ($4>=5013 && $4<=50130) print;}'
$grep -v '^@' tmp.sam | awk '{if ($4>=5013 && $4<=50130) print;}'| wc -l
17645
$grep -v '^@' tmp.sam | awk '{if ($4>=5013 && $4<=50130) print;}'|\
head -10 | less -SN
- 把sam文件按照染色体以及起始坐标排序
$grep -v '^@' tmp.sam | sort -k3,4 -V
$grep -v '^@' tmp.sam | sort -k3,4 -V |head -20 |less -SN
$grep -v '^@' tmp.sam | sort -k3,4 -V |tail -20 |less -SN
head -20结果
tail - 20结果
- 找到 102M3D11M 的比对情况，计算其reads片段长度。
$grep -v '^@' tmp.sam | awk '{if($6=="102M3D11M") print length($10);}'
113
- 安装samtools软件后使用samtools软件的各个功能尝试把上述题目重新做一遍。
另开一篇文章写作业。
